{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison to iPRG2012 consensus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "\n",
    "src_dir = os.path.abspath('../src')\n",
    "if src_dir not in sys.path:\n",
    "    sys.path.append(src_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import math\n",
    "\n",
    "import Levenshtein\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import squarify\n",
    "from matplotlib_venn import venn2, venn2_circles\n",
    "\n",
    "from ann_solo import reader"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot styling\n",
    "plt.style.use(['seaborn-white', 'seaborn-paper'])\n",
    "plt.rc('font', family='serif')\n",
    "sns.set_palette('Set1')\n",
    "sns.set_context('paper', font_scale=1.)    # two-column figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "psms_consensus = pd.read_csv(\n",
    "    '../../data/external/iprg2012/iprg2012ConsensusSpectrumIDcomparison.tsv',\n",
    "    sep='\\t', header=0, skipfooter=4, engine='python').rename(\n",
    "    columns={'bestSequence': 'sequence_consensus',\n",
    "             'Precursor_z': 'charge_consensus'})\n",
    "# get the same PSM identifiers\n",
    "psms_consensus = psms_consensus.set_index(psms_consensus['Index1_de'] - 1)\n",
    "    \n",
    "psms_annsolo = reader.read_mztab_ssms(\n",
    "    '../../data/processed/iprg2012/brute_force/bf_oms_shifted.mztab')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "psms_annsolo['mass_diff'] = (\n",
    "    (psms_annsolo['exp_mass_to_charge'] - psms_annsolo['calc_mass_to_charge'])\n",
    "    * psms_annsolo['charge'])\n",
    "\n",
    "psms = psms_annsolo[['sequence', 'charge', 'search_engine_score[1]', 'mass_diff']].join(\n",
    "    psms_consensus[['sequence_consensus', 'charge_consensus']], how='outer')\n",
    "\n",
    "# don't disambiguate between I and L\n",
    "for label in ['sequence', 'sequence_consensus']:\n",
    "    psms[label] = psms[label].str.replace('I', 'L')\n",
    "# remove SpectraST modification masses\n",
    "psms['sequence'] = psms['sequence'].str.replace(r'n?\\[\\d+\\]', '')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def edit_distance(seq1, seq2, normed=False):\n",
    "    if not pd.isnull(seq1) and not pd.isnull(seq2):\n",
    "        dist = Levenshtein.distance(seq1, seq2)\n",
    "        if normed:\n",
    "            dist /= max(len(seq1), len(seq2))\n",
    "        return dist\n",
    "    else:\n",
    "        return math.inf\n",
    "\n",
    "psms['edit_dist'] = psms.apply(\n",
    "    lambda psm: edit_distance(psm['sequence'], psm['sequence_consensus']),\n",
    "    axis=1)\n",
    "psms['edit_dist_norm'] = psms.apply(\n",
    "    lambda psm: edit_distance(psm['sequence'], psm['sequence_consensus'], True),\n",
    "    axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get unique keys for spectrum ID - peptide assignment combinations\n",
    "set_consensus = set(psms.loc[\n",
    "    psms['sequence_consensus'].notnull(), 'sequence_consensus'].items())\n",
    "set_annsolo = set(psms.loc[psms['sequence'].notnull(), 'sequence'].items())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASsAAAD5CAYAAAB/JRMkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd0VMXDxvHv9vTeCb2X0EMRRJAiCCqgiAUULCj2Cgg2REQQVFBERf0pqEgRQToISK+hI53Q0nvbvjvvH0heQEpQkns3mc85HE2y2X0SyJOZuffO1QghBJIkSSqnVTqAJElSSciykiTJI8iykiTJI8iykiTJI8iykiTJI8iykiTJI8iykiTJI8iykiTJI8iykiTJI8iykiTJI8iykiTJI8iykiTJI8iykiTJI8iykiTJI8iykiTJI8iykiTJI8iykiTJI8iykiTJI8iykiTJI8iykiTJI8iykiTJI8iykiTJI8iykiTJI8iykiTJI+iVDiBVbEVFRSQnJ5OSkkJycjKpqamYzWacTic2uwO73Y7T5USn1WE0GDCZjOj1eoxGIxEREcTExBAdHU1MTAzBwcFoNBqlvySplMiykkqV3W7n4MGD7Ny5k8OHj3DyzFnOnksiNTWVrPRUnE4nvoFheAeGYvQPQe8bgltnxCm0CI0Ojfb8H4Qb4XYi3G40woUeF8KSh6MwG2t+FubcTJwOG8FhEURERlE5thLVq8RSq2YNmjdvTvPmzQkICFD62yH9Bxp5+3jpZrm4mNZv3saOHTs5efwIAWExeEfVQhsUi94/BKNfKAb/EIz+oehMvjdtNORyWHEUZOMozMZemI2jIBtNYRqW1OPkJB0nLCKK5s1b0OGW1sTHt5QF5mFkWUn/ms1mY+3atfy+aDHrNmzi+LEj+IZG4xVRE2NETXyja+MTVROd0UvpqAi3C0vmWYqSj2JNO4Ez8yS5SccJj4ymTat47u51J3feeSdhYWFKR5WuQpaVdEMyMzNZtGgxs+b+xsb1awmKqYFXtZZ4V2qgmmIqqQsFVph0GOfZ3WQcTaBug4b0v7cPffv0pm7dukpHlC4iy0q6riNHjrBgwUJmzZvP4b8OElSjGb41WxFUKx6Db5DS8W4at9NOfuJe8k9so+DYdgL8/ejb5x769e3NLbfcgl4vl3iVJMtKuqKsrCy+/fY7pn75Nbn5+QTUboNvjVYEVGuMVm9UOl6pE0JgTj1O3rFt2BJ3YM1L56EHH+T5556lQYMGSserkGRZScWEEOzYsYOJn37G4kW/E1K3DYGNu+MXW7/CnxJgy00je99KsvasoF7dugx/9UV69+6NwWBQOlqFIctKwmw28/PPs/jo0ymkpmcR1KQ7YU27YfAJVDqa6rhdDnIObyZv71KcuSk8/dRTPDv0KSpVqqR0tHJPllUFdvz4caZ8PpX/ff8DPjH1CG56J4G1WqDRyAsbSsKcfors3UvIPriO2zp24vWXX6BTp04VfhRaWmRZVUBnz57ljVFvsfD33wlt0o2gJj3wCo5SOpbHctnMZO5fQ8HeJVSrFMmnkz6iffv2Sscqd2RZVSBZWVm8Pfp9vv/+f4Q160F4m/vQe/kpHavcEG4XWfvXkLHpJ1o0a8rkSROIi4tTOla5IcvKAwghyMrKIjk5maysLBwOB06ns/iPTqdDr9djMBjQ6/UEBwcTHR1NeHg4Op2OoqIiJkycxKSPPyGgbnui2j+I0T9U6S+r3HI77aQnLCFjy1y6devGJx+No1q1akrH8niyrFTA7XZz9OhRdu3aRWJiIkmnT5N86hTJSUmkpqeTlpODj8FApJcXIXoDRq0WvVaDXqNFqwG3EDjdAqdw4xSCXKeTVKuVPJsNH5MJs92Bj18gwdXicAfFYgiMwDe6Ft7hVc9fdyeVCqetiLSt88natZiBAwbw/uh3CA8PVzqWx5JlVcYuFFNCQgI7t25l58aN7D18mBCDgTiTF9XcbiKcLiJ0WqJ0OiK0OsJ1OrxvcNF2tdXCW3YHYSFR3FmrFb5GbzLNeWSZ80mx5HE44wypuZlERFdHF1UbTUQt/GJqywIrBY7CHNK3zCb3r3WMeP01hg8fJk8w/RdkWZWBo0ePsnDBApb+Op+EfXsJMZmIM5poZLPRWG8gzmgg+CYVRJ7bzbtOF5udLl5r35/WVepf9bFFditHM89xJOMMR3OTOZRxltTcTMJjaqCt2oKAOm3wiaopj27dJNacFNJWfUGw3s6cWT/SqFEjpSN5FFlWpcDpdLJlyxYWzJ7Not9+oyA3ly56A130eloYjTetmC632mphhN3BLTWbMbTV3fgYTDf8HEV2K4fST7Ml6TDrT+3H7HIRVLcNxpqtK8zZ66VJCEHG7uWkrZ/Ba6+8zFuj3pCjrBKSZXWT2O12lixZwm+zZrFk6VJi9Hq6uAXdTF7EGQylOjrJc7t5Vwg2252M6PAgLSrVuSnPK4TgdE4aG07vZ2PSIU5kJBFduzmiRmtC6rX3qIuW1caWl07Sss8INdqZK0dZJSLL6j86e/YsX06ezLfffEMNg4GeDiddTV5UKqPflqttVobbHLSr0ZShrf/daKqkss0FbDlzkHXnDrL73HFCG3cmqHlPvMMql9prlmdCCDL3LCdt3Qxee/Vl3hwpR1nXIsvqX3C73axevZrPJ0xg/caN9DGaGGjyok4ZXieW73bzjs3GFrQMv4mjqZJKLcjm9yNbWHhoC/5R1TE17UlQnTZodfKH7UbZ8tJJXv4ZYUYHc2b9SMOGDZWOpEqyrG5ATk4O33//PV98/DGmwkIGoqGPlze+2rK9POWk08Fgs4XGtZrxbHzpjqaux+Fy8ufJPfx2dAunczMJbn4ngU27y/O4bpAQgqw9y0nfMJPvv51O3759lY6kOrKsSqCoqIhPJk3ik4kT6ejjyyMuFy0MRkWOkq2zWnnBauXJ1ndzd/1byvz1r+VEVhILjmxhxdEEQpp0Jaxdf3kx9A0qSj7Gmd/G8vzQIbw3+h20ZfyLUM1kWV2Dw+Hg66lTef+dd2mt1fKa0UR1hdYUhBBMdwu+KCrivS6DaRpTS5EcJZFlzueHvX+w/OhOwlv3JrRVH3RGb6VjeQx7YTZn5o8lvlEt5v3yE76+vkpHUgVZVlfgdruZPXs2b77yClUsFkYYTTQyKHfI3iYEI1xu9mr0jOv2JNH+IYpluRHn8jL4ZtdytiUdI6L9AwQ164FWJ/d/Kgm3087ZZZ/ja05i9YolVK1aVelIipNldZkVK1Yw4uWX0aamMlyro51J2cPz6S4XTzgcBIVVZdRtD+Gt4PrUv3U08xxf7VzCsfxMAm99hNBGt8ltaEpACEH69gXk7pzP77/9yq233qp0JEXJsvpbRkYGzzz5JLvXr2e4VkcPo0nxM7f32e08YbHQM64jjzbritbDf8B3JR1lasJSso0+RPV8Ga/gaKUjeYTcEwkkLZ7ExPEf8NRTTykdRzGyrIA5v/zC8089xb16A6+YvG74OrzSsMlm5Wmzhddve4iONZooHeemcbndzDm4ju93rSK8wwDCWvaSo6wSsGSd49z8MTz75CBGv/uO4r9IlVChyyojI4Ohjz7K3nXrmeTjQwujOqZYa61WXrRaGdP1cZpVqq10nFJxOieNDzb8Qq7Jm5A75SirJByFOZya/SaDHuzLxAnjK1xhVdiymjt3Ls89+ST3arS84u2jitEUwAqLhddtNsbdMYS46BpKxylVLreb2fvWMmPvaiJuG0BICznKuh6HOY9Ts9/igbvv4PPPJleowqpwZZWbm8uTgwez588/maQ3qGY0BbDUYuENu50JPZ6mfkTFOfpzKieVDzb8QqbRh6i7XscUKPd8uhanpYAz896ld7db+WraFxWmsCpUWR05coS7ut1Be5uNkdob3yOqNP1htfCqzc7EO5+hbnjFu9bO5Xbz04E1/LxvA7H3jsK/srw337U4bUWcmfMW9/fqwmeTP60QhVVhymrpkiU8+sADDDeaeNBLXScorrdaec5qZXyPp2kYWU3pOIrafPogY9b9TPjtjxHa9A6l46ia01LAqV9GMfiBPkwYP67cF1a5LyshBJPGjmXi2A/40tePeJN6pn0AW2w2hpjNfNB9CE2iayodRxVO5aQyfOW3eNdtQ9DtT8idS6/BYc4j8ec3eGHIo7z7zttKxylV5bqsrFYrTwwYwL4VK/jGx6/Mtm0pqZNOB30KCnin6+O0jK2rdBxVybeaeWfND6QaTYTd8wZ6b3+lI6mWvTCbUz8N5+MP32PQoEFKxyk15baskpOT6X3HHUSdO8fHXj54q+yC0Hy3m15mM/1b9uSuem2VjqNKTreLadsXseLMX8T2exfv8CpKR1Itc8ZpEn9+g1XLl9CmTRul45SKcllWx48fp3P79vR3OHnB5KW6ubxLCAY5XYRWqscrt9yndBzVW3p0O59u+52q/d7FL7ae0nFUK+foNjJXf8HeXTuJjY1VOs5NV+7K6vDhw3Ru354XBAxQ2UL6BWOEIMHgz8c9hqLXyfWYkth0+iDvrfuZ2HvfIqCq3AL4alI2zcE7NYHdO7bg7a3Of///lrrmRv/R/v376dS2La+7hWqLap7dzhKrnTFdBsmiugHtqjbk/U4DSfp1DHkndysdR7WibulHji6EBwcOopyNQ8pPWR08eJCut97Km2jp5+2jdJwr2m23857ZzId3DCFQ3rb9hsVXrse4LoNJWfAheYl7lI6jShqNhqo9X2LD9r2MHfeh0nFuqnJRVkePHqVrhw6M0uro7aPOokp1uXjSYmH4bQ9RI0ReB/dvNYupxQedB5Hy2zgKzhxQOo4qaQ0mqtz7JhMmfcrixYuVjnPTePya1cmTJ7mtbVtecrl5QOG9p67GKQR9zWZaNezAI027Kh2nXNh25hDv/Pkjlfu/h18ledrHlRScO0TSb++TsGMbNWt6/jl8Hj2yys7OplunTgxFo9qiAvjCbscYEsPAJl2UjlJutK5Sn7due4Czc97Fmp2sdBxV8o+tT2jr+3hwwKO43W6l4/xnHltWTqeTfnffTReLlUf16t0q97DDwXSLhWHt+6vuFApP165qHEOadyNl3mictiKl46hSWMu7OZmWzyeTpygd5T/z2LJ65dlnEQf/4g0VF5VTCF52uRjS6i6iPGTfdE/Tp0F72kVUJXvRRwjh+aOHm02j1VHpzhd5593RHD9+XOk4/4lHltX0adNY+tNPfO7ljV7Fo5UvnE58AiJUd8us8kSj0fBS274EFWRTsGGm0nFUyTs0lrA29/OAh08HPa6sNm3axMhXXuE7Xz8CVXYJzcUOOxxMLypi+K0PyOlfKTPo9LzfZTC5e/8g+8BapeOoUmSrezidUcCnHjwd9KijgWfOnKF1kyZM0Oq4XaUnfQI4hOBuu50747pwT4N2SsepMI5nJvHc0qlUeWAsvjHlczvo/8KSdY4TM19jT8IOatVS730nr0a9Q5PLWCwW7unRgyd1elUXFcA0mw0ffzn9K2u1wioxon1/zs57D0dRrtJxVKd4OvjwIx45HfSYsnpr1CgqZ2YxREXbEF/JCaeD6VYLw2+VR/+U0LFGE+6p25LU5Z+Vu8tNbobIVvdwNquIadOmKR3lhnlEWW3ZsoUfp0/nfY1W9QXwoVvwULNu8uifgh5r2g1TxmlyD/6pdBTV0Wh1RHR+mjffGU1BQYHScW6I6svKYrEwqH9/RhuMhKr8wt9ddhu77A7ua1ix75yrNKPOwKhbHyRt1Vc4CnOUjqM6vlE18K7ShA8nTFQ6yg1RfVm9+dpr1M3Pp5dKL06+QAjBB1o9g5p3x6Q3Kh2nwmsQWZW7a7cid9Xncjp4BRHtH2by5ClkZGQoHaXEVF1WW7Zs4cfv/sf7XuouKoA/bTZSXYKe9crnLo2e6PH4Hoi002TL6eA/eAVHE9jgNka+PVrpKCWm2rKyWCwMur8/73l7q3765xaCcWh4smVP9PLmBqpxfjr4AOlyOnhFUe0e4Kcff+TUqVNKRykR1ZbV22++ST2rVfXTP4DfrRY0Xn50rNFU6SjSZRpEVqN3/bakrPhc6SiqY/ALJrR5T15/Y5TSUUpElWV14sQJvvvqK97TqetuNFdiF4KPnC6eatlL9UcqK6rBTbtB6nHyT+1TOorqRLa9l+XLlrN//36lo1yXKsvqzWHDeMzbhzCVT/8AZpmLiAmvLG+lpWImvYEhzbqTt+47udh+Gb3Jl/C2/Xj19eFKR7ku1ZXVnj17WLtiBU+qeDeFC1xC8IXLxePNuisdRbqOrrVb4GM1k3dks9JRVCeseU+27UjgwAF177yqurIa/tzzPGcw4qvii5QvWGuzEhwYToMKfst3T6DVaHm6ZU/y1v+AcLuUjqMqWr2RoMbd+Hiyutf1VNUIf/75J4d37+JhD1hUB/jeYKR3XXn9n6doW6UB0SZvMveuUjqK6oQ268Hs2b+o+qx21ZSVEILXhw7lNaMXRg9YqD7ldLK3oIDOtZorHUUqIY1Gw9MtepKz8SfcDpvScVTF6B9KQPWm/DBjhtJRrko1ZbVgwQIsKSn09pAbM/4oBD3qtZFnq3uYuKgaxEVUIW3nIqWjqE5QkzuZNFm9Z/yroqyEELw7fATDdHq0HjCqsgjB7MJCeteXe1V5oiHNupO9ZR5up13pKKoSUK0xuYVWNm7cqHSUK1JFWW3ZsoWitFQ6qfgONRdbbLVQP7o6sYHhSkeR/oUaIdE0iqpC1l/rlY6iKhqNhsAmPfjoE3XuJqqKsvp8/AQGarQeMaoC+EFnoHc9OaryZH1qt8G5Z4nSMVQnrHEX/li5ktTUVKWj/IPiZZWens7SlSu4z0OOAO6320lzOGhbpaHSUaT/4JaqjbDlZlCUckzpKKqi9/IltGEHvp4+Xeko/6B4WX375ZfcaTIR5AHnVQH85nDQo15bdB6SV7oynVZL73ptMe9ZqnQU1Qls1IXvZ/yodIx/UPQnzuVy8eWUKTxiUPdWxRf7A7i1SiOlY0g3wV312pJ2YANOa6HSUVTFt1IdMrOyVXefQUXLatmyZYQ6HMQZPePw/wmngyIBdcJilY4i3QQhPgG0q9FIniR6GY1GS0Dt1ixY+LvSUS6haFlNnTiRRzxo/6eVbsEtVePk7grlSJ86bSnctUTezfkyvtXjmTV3vtIxLqFYWaWnp7Nl23aP2K/qghV6A+3lwnq50jiqBj4aQVHSUaWjqEpAjaYc3LeHnBz1bFqoWFktWbKEDoEBeHvIKCXL5eJwTjbNK9VROop0E2k0GjpWboT1xDalo6iKzuBFUPUmLFu2TOkoxRQrq4WzZ9PZ5jlnEK+xWYmvWh+TB2xdI92Y9lUaYj0uy+pyPjVb8fOcX5WOUUyRsrJaraxdt47OXp5xxjrASoORdrENlI4hlYIGkdUoysvClqu+EyGVFFy7NX+u+QO7XR2DCkXKas2aNTT08SXYQxbXrUKwMT+PtlXlelV5pNNqaV+tEblHtyodRVUMfsEERFZl/Xp1XJakSFkt/OUXOjudSrz0v7LHbqdaWAzB3v5KR5FKSbvYBoiT25WOoTqmKs1YuVIdp3aUeVm53W4WLVpENw+aAu51OakbVkXpGFIpiq9cj/Qzh3Bai5SOoiqmyFqs36KO9bwyL6tdu3bh6xbU8KCF6v3evtQJqaR0DKkU+RhMNKtSj7wTCUpHURXf6Noc2LdXFXtclXlZrV69mo4ecsb6BfssFuqFy5FVedc6shb207uUjqEqBr9gtHqjKm6EWuZltWPDBho7PWfD/kK3m5SiAqqFRCkdRSpl9cOr4E5V1/VwauAXU5udO3cqHaPsyyph504aGzxnCnjA4aB2ZGV5W/gKoHZYJTLSzsodRC+jDa3Blm3KH3wo07LKzs4mKzeXGnr132n5gr1uF3VC5IXLFYFJbyQ2NBJz+imlo6iKb3QtNm3doXSMsi2rXbt20dDP32N2BAXY5+1LnVBZVhVFvbDKckO+y/jG1Obgvj2KL7KXaVnt3LGDOA86vwpgv7lILq5XIHWDY9Cky3Wrixn9QtAajCQmJiqao0zLaseffxKn8ZwdNs1uN8mFcnG9IqkbXhmbHFn9Q3DluuzevVvRDGU7Ddy9x6MW15NdLiIDg+XiegVSOzSWTLnI/g9avwjOnDmjbIayeqH8/HzSc3M8anE9ze0izC9Y6RhSGfIyGIkJicCSeVbpKKriMgVy9lySohnKrKySkpKI9vXzqMX1dJebUJ9ApWNIZSzaPxh7fqbSMVTF4B9K4pkKUlYpKSlEGjxnVAWQqtUQbPRTOoZUxkK9/HEUZisdQ1WMfiGcOXdO0QxlW1Z4zqgKIN3Hj3CfAKVjSGUs1OSPS5bVJQz+IaQpfOPTMiur5ORkwj3stIU0INRXTgMrmjAffwxmWVYXM/qHkp2RpmiGsiurM2eIdHnONYEAaQ4HoXJkVeGE+QYhimRZXUzn5YfDbsNisSiWoewW2BMTifCwUwDSrVbC5Miqwgn1CcBeIMvqYhqNBr/gMFJSUhTLUHZrVueSiNR5zgmhAOnmQsLk0cAKJ8w3kML8LKVjqI53YBjJycmKvX7ZlVVaKpE6zxlZ2YTA7nLha/ScHU2lmyPUJ4CCglx549PLGH0DycjIUOz1y+xcgkKLBX+vkt/QNFWjYarJgLeAkX/fsmuXTstSvZ7qbjfHdFpetdrxBz43GvAFdEKQpNXyss2OF+AG5hr06IFErZZwt5tHHSVb5HcIgUGnK5d3Xy60W3hr28/0rtGa2yo1wuZy8HvidramHmVS+8HFj8uyFvDzkfVUD4zkVH46vaq1pFpABIeyz/Hz0fX46E0A3B4bR+uo8nM/RYNOj16rw+20ozOUzi8r4bRTtOUH9GHVcBVmofUJwbthN9xFOdhObEZjMOFIPYJ3097ogytRtGM27vz/X+D2u/VJNEZvLPuX4rzoWkaf+P7oAiJLJbNbo8PhcJTKc5dEmZWVw+W6oRc7rNPSyulm/0VTx6+NBl62OajrdvO5xsBKg557HU4ChWDg3yX0qcnAMoOePg4nK/Q6QoSgq9OFABK1JS8eJ6DXedZ5YSU1/+RWqvpHFL99PDeFZuE12Jxy+JLHrTqzh1pB0fSo2pyE9BP8nridF5r0AuChOh2oX463ztFqtQh3aR4QEhgqxWGq0QYh3OT9NhJTrXaYE+bi2/ZRNAYTxqotQHd+V12tlz++8f2v+Ez+nV8oxZwXJdbocCp4RL/MfhpdLhe6GxildHS6WK6/dNoYLCD/76co0mio/vfRxYEXjZYE4PP3VhZr9HpauVz8atCTr9HQx17y3wouIdBpPWuNrSTWnNtHfEQtNl5UTA1Dq5Bhyf/HYwOMPhQ6zh/9KXRYqeofXvyxzamHOVWQjt3loHPlJviV0ghEKTqdDtylNw3U6E2YarQBQFgL0Rh9wGnDbcnFlrgNXDY0Bh9Mtdqdf4zThuXgCkCD1icIU/VWxc9lObgSjVYLOgOmWree///SoFW2rMrsp/FGR1ZX8oTdzi8GA18aDeRroKbr0n9MGRoNqRott/+9bXKGVoMNuNfhpKnLxTgvU8nzAvpyVlZnCzPJt1uoH1K5RI/vFBvH6fwM5h3fzOaUwzQNrwFAJb8Q7qneih5Vm1MjMIrpB1eWZmxF6LU6hKv0fzBtidso2vw/vJv2xm3Nx5WbjKFSI7zqd8WZmYj99PkbWBirtcSrfhe8G3bDkfIX9rPnd0AwVmmGV92OeNXvgrDkYTu6ttSyahSeBpbZT6MQ4j+dv+4ARnmZeNNm42m7g5YuN1+a/v/GE7nAdKOBkVYbF/Z18BGCun//dmzounRKed28AOVsvWp3RiIaNCw+tZNzhZnsyzzN9rSjV3389IMruT02jvtq3cLAeh35bO8S4PyIK8z7/PlndYMqcThH2WvGSoNGoymTBXZT9db4dXoO68HlIAQakx863xAA9GHVcWaePP//wZXR/H3qjyGiVvE6lS4wGo3+/M+B/qL3lwqNBncpjjavp8zKSq/T8V9WAOycn/oF/L1ZYYhbcKHjszUwzWTkOZudQGDT30cdm7rcpP69f1a6RkMld8l3OtSDon8xpeHu6vHcU6MVvaq1JNYvjMZhVWkVefWF8WxrIYEmXwACjT443OdHGotP7aTIYQUgw5JHhHf5O73D5XahKcU1S1deCs6s0wBotDq03oG4zblodAbE31NvtzkHrV8YAOY9C///cwszr/v+UuF2YVBwi6cyW7PSa3XcyKB6m07LFr2OZI2W3/V67nY6eczu4GOTkRi3m5NaLYP/XoMa8fcNU9/xPj/Nq+Ny087l4iG7gy9NRvI1es5ptbxus5U8rwacHnbGfUltSjnEucJMLE4bUT7B+Bm82JD8F1ang99PbueOqs0w6Qz0r92e5ad3UcU/nFRzLoPrdwYgzCuAH4+sI8Y3hOSibB5v0EXhr+jmc7ndxSOZUqHVYzu6CmdgNG57ERqvQAyxjdGY/DDvXojWJwi3rRCfhncAIKwFmPf+jkarR9gtmOJ6nn8ejQbzrl/RmPxw56fh3fzeUoss3C70Cm7xpBFltLFyRGAgK719CfeQc63y3W7is7NZ9cREpaNICuj0zavEvfILOqO30lFUI2nReD4aNoQHHnhAkdcvs2mgUa/HpoK7upaUgfNTAaniEULgdDpLdRroiXQIRaeBZVZWEaGhZHrQGpBJo8HlduMogyNCkrrk24rw8vJBq/OcLbjLgttuxs9Puf3dyqysoqOjSfOgNSCtRkOotw/Z5gKlo0hlLLMoD/+AEKVjqI4lL4uYmBjFXr/MyqpS1aqke9i0KsLbh0xzntIxpDKWWZSPl78sq8sV5WYQHR2t2OuX3ciqenVSPWjNCiDSaCSzSJZVRZNpzkPnH6p0DFVxO+04bBZCQ5X7vpTdyKpSJTKMxus/UEUidDo5sqqAsorycPnIuxpdzF6QTWBIuKIX9pfpmlX6DVxIrAaRdhtZRf+8Zk4q3zKtBbh95DTwYo7CbMIjlb3Zb5mVVUxMjEctsANEWi1kWeUCe0WTbSvEKKeBl7AXZFFJwcV1KOORVarVWlYvd1NE6nRk22RZVTTp5nwMfnJkdTFHYTZVK1fVuYFbAAAgAElEQVRSNEOZlVVUVBQ2l4scDzoiGKnVkVWYq3QMqQwJITiVkYR3qLI/mGojzDnUqFqy3TpKS5mVlVarpUm9euy7gT2llBap05GWL28cUJEk5WdiMPlg8JML7BfTWnMVPW0ByrCsAFq2b89+p+eUVfjfu0XK0xcqjiMZZwmqVFvpGKpTlJZIw4YNFc1QtmXVti37Fby26EZpNBrigkI4knFG6ShSGTmcdRZXeE2lY6iK22knN/U0jRs3VjRHmZZVixYt2O+wl+VL/mdxLieHM84qHUMqI0eyk/GOliOri5nTEqlctQbe3sruQFGmZVW7dm1y7HaPWmRv7HRwLEe5e6VJZUcIweHUU/jKsrpEUcoxWrRsoXSMsi0rrVZLs4YNPWqRvbHByOG0U0rHkMpAUn4mRpO3XFy/jC3tBLfd0lrpGGVbVgAt2rdnvwfdPLKyTofdYZeL7BXAkYyzBFYqP/c/vFnsGSdo2bKl0jHKvqxuad+ebXrP2C0Uzi+yNwoKlovsFcDezFM4o+opHUNV3E47+WlnFF9cBwXKqlu3buzIy6PQgzbiayIX2cs9IQTrT+4juE4bpaOoiloW10GBsgoICKBNs2ass3nOpTdxTgfH8uQie3l2LCsJt8GIV5iyZ2mrjTn1OG1bxysdA1CgrADuGTCAP5R44X+phdHEvrNHcXrQUUzpxmw8dQC/2m0U3QJFjdypf9GxQ3ulYwAKldXdd9/NGqsFl4dsxhet0xHr7cP+lJNKR5FKyaakv/CupfwRLzVxu5ykH9lBz549lY4CKFRWVapUIbZSLAl2zzlBtJsQbDp7UOkYUinIKMzlbE4GfpWVvZxEbQrOHKB6jVqKXxN4gSJlBXBX//tZhWeMrAC66nRsPLWfMrrNolSGNp0+QEy9VmjlrbcukX9sK/ff21vpGMUUK6t7evfmDw9aA2pkMGC3Wzmdm6Z0FOkm25h8GHf1VkrHUBUhBIUnttO3jywrmjdvjs1oZJ+HTAU1Gg1dAwLZeGq/0lGkmyjbXMCec8cJrKn8SY9qYsk4jUmvpVGjRkpHKaZYWWm1Wp567jlmOmxKRbhh3axmNp/7S+kY0k20+Og2Qhq0R+/lq3QUVck9upXe99ytqqOjipUVwBNDhrDUaiXXQ04QvcVk4njaWXIscqvj8sDldvPboU0ENVfH0S41sZ/ayf339lE6xiUULauIiAh6dO3KPHORkjFKzKTRcKu/P1vOyKOC5cGWMwcxBITKXRYuYy/MpjDjDB06dFA6yiUULSuA50aMYIbbhdtDjrL1cLlYe3qv0jGkm2D+4c2YmspR1eVyD22kR487MarsPp+Kl9Utt9yCb1g4m+yesXbVw9ubQ0knOJeXoXQU6T84l5fBX+lnCGl4m9JRVEUIQe7eZTz/7FClo/yD4mWl0Wh4ZtjrzPCQdStvjYZ+/oEsPLRZ6SjSf7Dg8GYimt2BVq+u0YPSCk7vI9DHyK233qp0lH9QvKwABgwcyFa7jSSnU+koJTIQwdLDW7A5PeO0C+lSVoedxUe24du0h9JRVCd3z1JeefF5VR0FvEAVZeXn58fjQ4bwmYf88FfT62nsH8Ca47uVjiL9C/P+2oBf1cZ4BavjMhK1sBdkkZe4h0cfGah0lCtSRVkBvPHmmyy12TjpIbfqGuSwseDIJqVjSDco32rmxz2rCbvtUaWjqE7m7uXcf39/AgIClI5yRaopq5CQEF595RUmWC1KRymR201eZOVncThd7iDqSX7at5qQ+u3wlvtWXcLtcpKzdzmvvPic0lGuSjVlBfDSiBEkaLXs9YBLcHQaDQMNBhYclgvtniKjMJcFf20mqP1DSkdRndyjW6lVqyZxcXFKR7kqVZWVj48Pb40dy3iXZ0wFHzQY+PN4AvlWs9JRpBL4btcKYlp0xxQQrnQU1Snct4zhr7yodIxrUlVZATzx5JMk+fiwwQO2PQ7T6ejm7cPcv9YrHUW6jtM5aaw9tQ+/NvcrHUV1Cs4dwpWXQt++fZWOck2qKyuDwcDYSZP4EOERe0e9YjDw6941ZJvl9YJq9s3eFYS27oPe21/pKKoihCBrwwzGjhmtujPWL6e6sgLo168fhkqVmO0BZ7VX0evp7ePLjD2rlI4iXUVC0lF2JCcS2ko9ezOpRd6JBEyuQgYNGqR0lOtSZVlptVq+nTmTcRYzyS71nyj6kl7PyiNbSc7PVDqKdBmzw8YHG2YT1f05dAYvpeOoihBu0tZ9z8cTxqHXq3+XVFWWFUDjxo15/uWXGW4uUv10MEynY1BAIN8mLFc6inSZL7cvwlitMcF15M0gLpd1cB1RIf6qX6u6QLVlBTDynXfIjohgjkX9R9ue0mjYceYgxzOTlI4i/S0h6SirTx8grPMQpaOojtvlIH3Dj3wxeZIqL625ElWXlcFgYMb8+Xxgs6p+Ouiv1fKctzdfJyxROorE+enfuPW/ENn9ebmofgUZu5bRuEE9OnXqpHSUElN1WcHf08FXX2W4w6766eBAg4HE9FPsST6udJQK78sdi/Gu0YRAOf37B5fdQuaWOXz2yUdKR7khqi8rgJFvv01OZCRzVH6hs0mj4XWdji92/I7LQ7a8KY8Sko6y6tR+gm6X078rydg2j66dO9GsWTOlo9wQjygrg8HAD7/8wgcWC0cd6j67va+3D77mPObs/1PpKBVStjmfMetnyenfVRSlniR79zI+m/yJ0lFumEeUFZyfDk76/HMet5rJUfGoRavRMMlgYOauFZzOkfcYLEt2l4ORa77Hu3E3efTvCtwuJ0lLP+HjiROoVKmS0nFumEaofSHoMi8/+ywJM2cy0y8AvYqPYnxntzPP5M/Uu15Ep/WY3wkeSwjBhxtmsxdBVJ9RaDTye3655PU/UUOXxtpVyz3mCODFPO5v9KPJk/Fq2JAxRYVKR7mmQQYDPrZCZu9fq3SUCmHe/nXszE4iotdrsqiuoCj1JLl7lvDTD995ZFGBB5aVXq9nzrJlrPf3Y5aKz7/SajRM1Gr4cddKOR0sZTvOHuZ/e1cTce876IzeSsdRnQvTv08mfeSR078LPG4aeMGRI0doHx/PdKMX8SaT0nGu6ju7jXleAUztJaeDpeFcXgZPLfyUqL4jCajWROk4qpS84SdqaD13+neBx/701K1blxmzZzPUbuWsim80MchgxMeSL48OloICm5lhq74l5LZHZFFdRVHqSXJ3e/b07wKPLSuAHj168Oa4cTxos6j2DHetRsMknY4fdy2Xl+LcREV2C6+s+Bp3rVaEtZA3Kr0Sl8NK6vJPmPzxRI+e/l3g0WUF8OyzzzL0jTd40GIm3eVSOs4VVdXrGWPy4o2V08mxyH2v/iuLw8ZrK6dTFFOPqM5PKh1HlYQQJC+bTOd2rTxi+5eS8PiyAhg2fDiPvPgiD5kLyVJpYfX28aG3Dt5a/T0OlY4CPYHVYWfEqm/JD69OxB3PePzUprSkbfqFUPL5/rtvys33qFyUFcBbo0dz39Ch3G8uJE2lhTXcaCLEnMunm+crHcUjmR02Xls1nfTgGMJ7vCBPUbiK7CNbKDywkuVLfsfLq/zs4VVu/rY1Gg1jxo/n4Zdfpp+5UJVrWFqNhs/0Og6cO8j8AxuUjuNRCm0WXl7xFRmhVYns+TIarU7pSKpkTkskadkUli1eSExMjNJxbqpyU1YXvDV6NE+PHEk/cxGJKjxK6K/V8p3BwP8SlrAr6ajScTxCjqWAF1d+RU5UXaLvfEEW1VU4zHmcnj+GLz6bTHx8vNJxbrpyV1YAr48YwRsTJtC3qICNKrxLTnW9nqleXry7+nu5FfJ1HM9M4vGFn1BYrSUxdzwjp35X4XY5ObfwQx5/5CHV3v79v/LYk0JLYu3atTzQuzfPabQM9vFV3ULjtxYzPxp8+OLul/A1lp+1hZvlzxN7GL9xDuF3DCW0kedsElfWhBCk/vEldfxtLF+6GJ2ufI48y3VZAZw8eZK7b7+dxtnZjPULwKiiwhJC8IbZzOHgSCZ0HYKXQd23QiorbuHm+4QV/HZsO9H3vo1vTG2lI6la2safMSTvZMfWzQQGBiodp9SU+7ICKCgoYEDv3qTt3MlX3r6Eqeg3j0sIXrbZSA6MZFzXJzDpDUpHUpTZYWPs+lkcsRYRee9bGP1ClI6kammb5yBOrmfb5o1ERkYqHadUVYgFAH9/f35btYquzzxDL3MhBxzq2XFUp9HwsclEaEEGb6+p2OdgpRRk8/TiKRz3DabSw+NlUV1HxvYFOI+uYdP6P8t9UUEFGVldbO7cuTzzxBM8bfLiSb0BnUqmhQ4hGGK14gyvyrudHsGgU/993G4WIQTLj+5g8raFBN/Sn8jWfVS3vqg2mQmLKdqzkK2bNlC1alWl45SJCldWAImJiQx+6CEKDx9mksFELYM6pl42IRhqtWIOjWXM7YMqxJQwoyiXjzbO4ZilkPCer+AbXUvpSKqXvm0+toPL2bj+T6pXr650nDJTIcsKwO12M23qVN4eMYKhegNDfP1UMcpyCMHzNiuZgVF80OXxcrvoLoRg+ZHtfLb9dyLjexFwywNodeW/nP+rtE2/4Dyxns3r/6Ry5cpKxylTFbasLkhMTGTQvfdSdPwEk7y8VTHKcgrBq3Ybp/xCGdflCfxM5WtDuYyiXD7aPI9jRflyNFVCQggyNsxEk5TAxnVriY6OVjpSmavwZQUXjbLeGMnTPj48rjdgUniU5RaCty0W1ulNfNjtSSoHRSia52Zwud0sO7adqdsWEdSiFxHt5WiqJFwOKynLphCmzWfV8qWEh4crHUkRsqwukpiYyAvPPMPezZt5Raenr5e34lPDmRYLE+123rr9EVpVrq9oln9LCMGm0wf4KmEphd7+hHcdim+0PHeqJGx5GZyZP4Yu7Vryw/++LVcXJt8oWVZXsHHjRoa/8ALZJ04wTG+gq8lL0aNTW202hlqtPNSsG/fHdfSoI2V7U07w5c4lpDvtBHR4lMA6bTwqv5IKzhzk9IJxjBz+Om8Mf73Cf99kWV2FEILFixfzxvPP45OdzRtGL1opuNf7WaeTwU4XtaJq89qt/TCqfPp0IiuJr3Yt46/MZEI6DCAkrrO8APkGZOxZQfq6H5j100x69bxT6TiqIMvqOlwuFz/OnMnbw4ZTRwPPuQUtDEZFfssVud285BYkGX0Y2+VxQn0CyjzD9ZzISuLnv9ax8dRfhN3Sn7AWPdHqy+cRzdIg3C6SV3+DJmUfK5ctpl69ekpHUg1ZViVks9mYPn06n44fj7fZzAA39PHywreM71jjFoJPbTZmudy83/Vx6kcof0Kgw+Xkz5N7+O3IFk7lZRLQ7E7CW/dGb/JVOppHcZjzSV48gboxQfw2bw7BwcFKR1IVWVY3yO12s3r1aqZOnMj6DRvoYzQx0Giidhmf8rDUYuENq4U+jW9nYNMuipzxnlqQze+HN7Pw8FYCoqphbNqLoDpt0Fags+9vluwjW0hd+QVPDH6ECeM/RK+X38PLybL6D86ePctXU6bwzddfU0urY6BGQ2eTF95lNNpKcbkY7nJzTmdk5G0PUycsttRf0+lysTPpCAuObSPh7BEimnTGv2lPvMOrlPprl0dOSwHJq75Ek3WCWT/+QIcOHZSOpFqyrG4Cu93Ob7/9xleffsrO3btp6+NDFzd08fIiopR3eBBCMMduZ6zFQt9GtzGwWdebPsrKt5rZeuYgm5IOseXUQfzDK+MV143QuE7yDsj/QfaRLSSvmMrDD/bn44/G4+srp83XIsvqJsvJyWHZsmUsmPkjK9f9SU2jkS5uQVcvb+rp9aW2MJ/icjHc7eac9uaMss7lZbDp9AE2Jh3ir5REoms2QVRvTVDtVhj9Q29S6orJaSkg6e/R1C9yNFVisqxKkd1uZ8OGDSz4eRaLFv2O22ajpbcPjaw24vQ64gxG/G/ilFEIwRyblbFWK33jOjKwaclGWVaHnWNZ5ziScZYjOcnsTz9FjqWIkLptMNVsTUCNZugMFfdkxJsp+8gWUlZO5cH+/flkohxN3QhZVmVECMGhQ4fYvn07O7dtI2HjRvYfO0aUlxeNdToaOV3EGQxU1xsI12rR/4cRWIrLxXCHk9M6A0Na9qJ9tTg0Gg0ut5tcayHJ+Znniyk3mUMZZzmXlUZ4ZBVM0bUQEbXxjamNT1RNud/5TWTOOE3Wxplo85L4aeb3cjT1L8iyUpDT6eTw4cMkJCSwY/16ErZv51RSEpl5eQR7exNl8iJCpyXc7SbS7iBCqyVUq0WPBoMGdGjQacAlwIXAKcCBINftJtXpZLcQ7HQ6cKPBqDdSaDXj4+NHQHAkhqiaaCJq4xtdG++IqvJcqFJiy0snbcNPFJ3cyaiRI3jh+ecq9CUz/4UsKxVyuVykp6eTkpJCcnLy+f8mJZF88iQZ6ek4nE5cTidOlwun04lOp0Ov06HX69HrdISEhRFTrRoxlSsTGRnJnr37+OLLr9EGVyGq46P4RNZQ+kss9xzmPFI3zSb3wGqGPj2Ut0aNKNf7o5cFWVYVhM1mY8rnU3l/7Dh8qzUlov0AvIKjlI5V7rjsFtK3LyBzx0L69evHuPdHV8jtXEqDLKsKpqCggA8nfMSUKZ8RFteJoBb3YAqSpfVfuWxmsvf9QcbWOXS+vRMfffgBtWrJfbpuJllWFVR6ejrjPhzPt9/9D9/Y+gQ27kFgrRZyUf0GmdNPkb17KdkH/6RTx46MGf0OzZs3VzpWuSTLqoIzm838+NNPfPTJFNIzcwhq0p2wpndgUOFF0mrhdjnIObyFvL1LceYm89SQITw79CliY0v/CoKKTJaVBJw/tWL79u1M/GQKS5YsJrhuW4Kb3olvTJ0Kv4/SBbb8DDJ3LSNn30rq1a3L8FdfpHfv3hhUsBV2RSDLSvqHzMxMvvn2O6Z8/gUunRfetdoSULsN3uFVK1xx2QuzyTu2DcepnWSfOsCDDz3Eyy88R4MGDZSOVuHIspKuyu12s3btWn6dv4Bff1uAzSnwrRlPQK3W+FeNK5e7KwghsGScJufIVqyJ2zFnnqNL127cf29v7rnnHvz8/JSOWGHJspJKRAjBgQMHmDd/AbPn/cbpxOOE1I7Hu0Y8QTVbovf2Vzriv+Z2OSk4vZ/CE9spOL4No15Lr1538dD9fenQoQNGozxhVg1kWUn/SkpKCosXL+aXefPZvHEDAeGVMITXwBBeE9/o2vhEVkdrUG4b6KsRwo01K4milGNYU4/jyjxJbtIxatWuy/339qF373to1KhRhZvueoJyW1bjx49ny5YtLFiwoPh9b7/9NosWLeLuu+9m9OjRxe8fM2YMK1as4Pbbb+e9997jk08+4ddff6VBgwZ89dVXxf9wP/vsM2bPns3tt9/Oyy+//I+dHF0uFxMnTiQtLY2goCDy8vKwWCy8/fbbREVd+1ymadOm8dNPP1GpUiWGDBlC586db+J3o3RZrVb27dvHzp07WbdpGzsTEjh76gQBEZUxRdbE+HeBeUdULdMtZdwuJ7acFIpSjmFOOY4r8wR5SccJCg6hafPmdGjbmlat4mnevDkhISFllkv6d8plWdntdh577DF27drFL7/8QuPGjYs/ds8995CcnMzIkSPp06dP8fsHDhzIzJkzL3k7KSmJvn378txzz131cRdbsGAB8+fPZ8aMGcXve+mll3j88ceJi4u7bu6BAwfSsWNHHn/88Rv6etXIarWyf/9+du7cycYt29ixM4FTJ4+h1enxDQrH6BeCxicY4RWE0T8Eg38IRr8QDH4haI1eaLV60OrQ/P0H4Ua4XcV/3E4HjqIcHAVZ2AuzcRRk4yrKAUsOrqIcLHmZWM2FREbH0LRZc25pHU9rWUwerfytkAIrV66kT58++Pv7M2vWrEvKKiAggFdffZUXX3yR+vXrX3ND/qlTp/LQQw/RuHHjEl0lX1RURE5ODoWFhcULsWPHji1e87Db7Xz00Uc4HA6EEHh5efHqq6+WaE0kLS2N8ePHExYWRlZWFm3btuW+++677ucpxcvLi/j4eOLj4xk6dChwft0rLy/vkmseT585x+lzSSQnHyP5cAoZaalYLBacF65/dDpxu5xotTp0eh06vQGdTofRaCQsPIIqUdHExsZQNbYeVSrHEh0dTXR0NDExMYSFhaEr5c0PpTIkyqGnnnpKFBUVibVr14omTZqI/Pz84o8NGDBACCHEl19+KTp37izy8vIuef/lj1u0aJFo1aqVOHv27BUfd7GCggIxYMAA0bJlSzFs2DCxdOlSYTabiz/+2WefiREjRhS/PXz4cPH5559f8prffPPNFZ97wIABYv78+UIIIRwOh+jWrZvYvn379b8ZklROlLtrK06ePElkZCQ+Pj506NCBoKCgS9atLnjqqado2LAhr7/+OuIaM+FevXrRt29fXnjhBWw22yUfe/TRR+nevTvdu3fH5XLh5+fHzJkzmTlzJjVr1mT69Ol06dKFI0eOALBmzRri4+OLPz8+Pp41a9Zc92sqKipi+/btxZ+r1+tp2rRpiT5XksqLcldWs2fPprCwkLFjxzJu3DgiIiKYPXv2FR87btw4kpOT+eKLL675nK+99hpBQUGXLMoD/PDDDyxfvpzly5ej0+lITU0lOzubevXqMWTIEObPn098fDxz58696nP/l6NO8oiVVJGUq7Ky2WycOXOGSZMmMWrUKEaNGsVnn33GqVOn2Llz5z8e7+PjwxdffMHMmTPJyMi46vPqdDo+/vhjtm3bxr59+676uM2bNzNv3rx/vP/CFiGdOnVix44dxe/fuXMnHTt2vObXNHDgQHx9fWnVqlXx1+ByudizZ891P1eSypNyczTQarXy2muvkZmZyZgxY6hduzYACxcuZMKECURGRuLn50diYiJ33HEHb775ZvHnbty4kdGjR7Nq1SoApkyZwty5c+natSsvvfQSAQHnL+o9fPgwDz74ILt3775ihkOHDjF58mSCgoLw8fEhNzeX4OBghg8fjtFoxGazFS+wAxgMBoYNG4bRaOTrr79m5syZxMbGUr9+/eLnXLVqFRs2bCAlJYUJEyYQHh5OVlYWrVu35v777y+V76UkqVG5KStJksq3cjUNlCSp/JJlJUmSR5BlJUmSRyiXZ7BL5dMTTzzBE088QZs2bZg1axbff/89BoOB5s2bY7FYyMrK4s0336RGjRqMGzeOOXPm0KxZM2JjY7FYLKSlpfHWW28VH3wBWLFiBcuXL8fX1xetVsvp06dp0KABzzzzDP7+/syYMYPdu3dTqVIlTp8+TePGjXnyySeLP/+PP/7g119/JSoqiuzsbN55553iy3ncbjdz5szhk08+YcqUKbRu3br48673vNIVKHpKqiTdgHPnzgmbzVb89vDhw8Xo0aOL3/7qq69Ev379it/u1KmTWLZsWfHb06ZNEwMHDix+e8aMGWLQoEGXXGWQnZ0tOnbsKI4fPy6EEOLhhx8WhYWFQgghbDabaNeunVi3bp0QQoiMjAzRtm3b4qsgfvjhB/H6668XP9evv/4q1q1bJzp16iS2bt16yddyreeVrkyOrCSPsHDhQr744gueeuop+vbte8XH1K1bl2nTpl31ObKysigoKADOn+ry8ccf8+233+Lt/f87QQQHBzNx4kQiIiIA+PHHH4s/ZjQaiYyMLD4nb8mSJcTFxRWf2tKpUycmTJjAmDFjMJlMV815vee9XEpKCuPHjyciIoKCggKcTicfffQRAHPnzmXr1q2EhYWRmZnJsGHDiIyM5NVXX2XlypW8+uqrbNy4kbNnzzJq1Cg6dOiA2+3m/fffB86fQ3ju3Dk+/vhjvL29SUhI4OeffyYiIoKkpCQefvhhWrduzYcffsi8efMYPHgwBw4c4OTJk5f8XXz++eekpqbi5+fHqVOneOedd9izZw+ffvpp8eO++uorvvzySxYtWkRsbCx//PEHS5cuJTIyktOnT/Pwww/Trl27q37P5MhK8hjDhw8Xv/766yVvXzyyeu+998TQoUOL3+7UqZMYPHiwGDFihOjatavo1KmTWL9+vRBCiL1794o6depcct3o9SQlJYlOnToVf86YMWPEyJEjiz9usVhEnTp1RGJi4iWfd6WR1bWe93IPP/zwJV/3I488IoQQYtu2baJr167C4XAIIc6P5C587MLr/vTTT0IIIdasWSP69+8vhBDi4MGDonv37sWPmzFjhigoKBA5OTmiffv2Ijs7WwghREpKimjVqpUoKioSQpy/PnXChAlCCCGOHDki2rdvL4QQIjc3VzRr1qw4x5IlS4q/B5f/nXXq1Kn4Otu77rpL7NmzRwghxPHjx8Xy5cuv+j0SQo6sJA+XkJDA2LFjMZvNBAUF8eGHH17y8fvvv5/u3buzbds2FixYQPv27a/6XPv372fu3LkkJiZy3333cc899xR/zGKx8O677zJ58mT8/a+8K+qFy59u5DKo6z1vUVERO3bsuOTr+uGHHwBYu3YtzZs3R68//2McHx/PyJEjKSoqwtfXF6B4naxGjRrFI7cqVapgNBoZOHAg3bt3p1evXvj5+bF27VpsNtsll59Vr16d9PR0qlWrdtXn8/f3p0mTJvTr149evXrRs2fP6+7fBtC1a1defPFF7rrrLrp3784dd9xxzcfLspI8WosWLRg1atR1H9e6dWu+/vrr4s0X69Spg6+vL0eOHKFly5YAxMXFERcXx8CBA8nMzCz+3IKCAkaNGsWLL75Iw4YNi98fGxvL5s2bi99OT09Hr9eX6Af1Ws/7X11clibT+d1adTpd8QX7fn5+LFiwgF27dvH7778zdepUZs2aBZzf2ufi76fFYil+joufT6/XFz+fVqvlf//7H3/99ReLFy+md+/eTJ06lRYtWqDRaHC73cWfb7fbi///+eefp1+/fixbtoznn3+eAQMG8Nhjj13165KnLkgVxmOPPcbUqVNxuVx4eaXlGiEAAAKWSURBVHnxyiuvMGXKFMxmc/FjhBCX/EBlZGQwbNgwXn75ZRo2bEhKSgrffvstAD179mTfvn3k5eUBsHr1anr06HHJD/fVXOt5L3b5daEAL774IoWFhXTs2JHdu3fjcrmA86PMVq1a4ePjc83XPnDgAAsXLqRFixaMHj2aZs2acezYMZo1a4bNZiveJcRutzNo0CCcTuc1ny89PZ2vv/6aBg0aMGzYMO6++24OHDgAQHh4OGlpacD5tbesrKzizxs7dixRUVEMHjyYkSNHkpCQcM3X0b377rvvXvMRkqQCixYtYvHixSQlJVGjRg3WrVvH8uXLSU1NJTMz85LTAgAmTJjAzp07ycnJwd/fn+rVq1OlSpXiBelKlSrRrVs3vL29mTZtGhs2bGDDhg3MmTOHBg0acO+99+Lv78+jjz7KgQMHWLRoEd999x2zZs2ifv36tG7dGl9fX6pWrcqnn37Kjh07OHnyJG+//Xbxgv2hQ4eYNm0ae/bsITs7m9zc3OKNIK/1vJdr27YtP//8M7t27WLZsmV07NiRRo0aERsbi8Fg4IcffiAhIYFDhw4xevRo/Pz8+Prrr9m4cSOZmZm0a9eOCRMm8Ndff6HT6WjQoAHTp09n//79rF+/HqPRyCOPPIKvry/NmjVj6tSp7N+/nxUrVvDEE09QuXJl5syZw8qVK0lNTSU+Pp6vvvqKhIQEbDYbrVq1YsaMGezatYstW7aQlZXF008/jclkIjY2lu+++479+/eTm5vL8ePHOXfuHK1bt2bz5s0sXbqUXbt2kZCQwNChQwkPD7/qvwF5baAkSR5BTgMlSfIIsqwkSfIIsqwkSfIIsqwkSfIIsqwkSfIIsqwkSfIIsqwkSfIIsqwkSfII/wcpBSQKGrdo9AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x311.496 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "width = 7\n",
    "height = width / 1.618    # golden ratio\n",
    "fig, ax = plt.subplots(figsize=(width, height))\n",
    "\n",
    "v = venn2([set_annsolo, set_consensus],\n",
    "          set_labels=['ANN-SoLo', 'iPRG2012 consensus'],\n",
    "          set_colors=[next(ax._get_lines.prop_cycler)['color'],\n",
    "                      next(ax._get_lines.prop_cycler)['color']],\n",
    "          alpha=1., ax=ax)\n",
    "c = venn2_circles([set_annsolo, set_consensus], linewidth=1.0, ax=ax)\n",
    "\n",
    "# plt.savefig('iprg2012_consensus_venn.pdf', dpi=300, bbox_inches='tight')\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "psms_match = psms[psms['sequence_consensus'].notnull() &\n",
    "                  psms['sequence'].notnull() &\n",
    "                  (psms['sequence_consensus'] == psms['sequence'])]\n",
    "psms_different = psms[psms['sequence_consensus'].notnull() &\n",
    "                      psms['sequence'].notnull() &\n",
    "                      (psms['sequence_consensus'] != psms['sequence'])]\n",
    "psms_unique_consensus = (psms[psms['sequence_consensus'].notnull()]\n",
    "                         .drop(psms_match.index, errors='ignore')\n",
    "                         .drop(psms_different.index, errors='ignore'))\n",
    "psms_unique_annsolo = (psms[psms['sequence'].notnull()]\n",
    "                       .drop(psms_match.index, errors='ignore')\n",
    "                       .drop(psms_different.index, errors='ignore'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# identical PSMs: 4185\n",
      "# conflicting PSMs: 742\n",
      "# unique PSMs ANN-SoLo: 1084\n",
      "# unique PSMs iPRG2012 consensus: 2913\n"
     ]
    }
   ],
   "source": [
    "print(f'# identical PSMs: {len(psms_match)}\\n'\n",
    "      f'# conflicting PSMs: {len(psms_different)}\\n'\n",
    "      f'# unique PSMs ANN-SoLo: {len(psms_unique_annsolo)}\\n'\n",
    "      f'# unique PSMs iPRG2012 consensus: {len(psms_unique_consensus)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAAEeCAYAAAAZ5BURAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAH+1JREFUeJzt3XtUVWX+x/EPYKjcBCU10cRbRxS1i4rNeBlva9Sxm2bW6DB2UZwmynTUnPI3VlaKmRResholu2mWtya1Mly5UkHN+wqQhFHB26gcVEBR2L8/XJ0V6engwbMPsN+vtVirvT1776+7s/j4PPvZz+NjGIYhAAAsytfbBQAA4E0EIQDA0ghCAIClEYQAAEsjCAEAlkYQAgAsjSAEAFgaQQgAsDSCEABgaQQhAMDSqnUQ2mw2b5cAAKjmqnUQAgBQWQQhAMDSCEIAgKURhAAASyMIAQCWRhACACyNIAQAWBpBCACwNIIQAGBpBCEAwNJqebuA6i4vopnbx0bkHbmBlQAA3OGRINy7d6+mTp2qcePGqXfv3nr//fd18uRJ+fv7q7i4WM8995wkKT09XcuWLVOdOnXUuXNn9evXzxPlAADglEeCMC8vT1FRUY7ty5cva+LEiZKkuLg47d27Vx07dtT06dO1YMECBQUF6cEHH1TPnj3l7+/viZIAALgmjzwjHDhwYLntxx9/3PHfZWVlqlu3rkpKSnTq1CmFhITI19dXjRs3VlZWlifKAQDAKVOfEf73v/9VWFiY2rRpoxMnTqhu3bqOPwsICJDdbnd6bFJSkubOnWtGmQAACzEtCO12uxYuXKhp06ZJksLCwlRcXOz486KiIoWGhjo9Pj4+XvHx8eX2sR4hAKCyTHl9oqioSK+//romT56sixcvateuXfL391d4eLjOnj2rsrIyHT9+XG3atDGjHAAAHDzSIlyzZo0yMzNVWlqq8PBwffDBB/rxxx8VHx+vixcv6uGHH9Ydd9yhF154Qa+//rpq166tJ598koEyAADT+RiGYXi7CHfZbDZlZmZ6tQbeIwSA6o2ZZQAAlsbMMl5EaxIAvI8WIQDA0ghCAIClEYQAAEsjCAEAlkYQAgAsjSAEAFgaQQgAsDSCEABgaQQhAMDSCEIAgKURhAAASyMIAQCWRhACACyNIAQAWBpBCACwNIIQAGBpBCEAwNIIQgCApRGEAABLIwgBAJZGEAIALI0gBABYGkEIALA0ghAAYGkEIQDA0ghCAIClEYQAAEur5YmT7t27V1OnTtW4cePUu3dvlZSUKCEhQeHh4crPz9ekSZPk5+en9PR0LVu2THXq1FHnzp3Vr18/T5QDAIBTHmkR5uXlKSoqyrG9cuVKRUZGauzYsapbt65SUlIkSdOnT9f48eM1adIkzZ8/XyUlJZ4oBwAAp1wGYWJionbv3q0vv/xS3bt319tvv+3ypAMHDiy3nZaW5gjG6OhopaWlqaSkRKdOnVJISIh8fX3VuHFjZWVlufnXAADAPS6DsHbt2rr99tv1/vvv6z//+Y8uXLhw3Rex2+0KDAyUJAUEBCg/P1/5+fmqW7eu4zMBAQGy2+1Oz5GUlCSbzVbuBwCAynL5jPDs2bPauXOnmjRpotDQUPn7+1/3RUJDQ1VYWChJKioqUlhYmMLCwlRcXOz4TFFRkUJDQ52eIz4+XvHx8eX2EYYAgMpy2SKMiIjQSy+9pLi4OG3cuFHHjh277ovExMQoPT1dkrRv3z7FxMTI399f4eHhOnv2rMrKynT8+HG1adPm+v8GAABUgo9hGMa1/uDMmTOqX7/+Vfvtdvtvttwkac2aNVq8eLFat26t2NhY2Ww2zZw5U/Xr11dBQYEmT57sGDX6ySefqHbt2oqJibnuUaM2m02ZmZnXdcyNlhfRzCvXjcg74pXrAkBN4zQI//nPf+qee+65av8XX3yhV1991eOFVQRBCACoLKfPCNevX699+/YpLCxMv8xKd7pGAQCoqpwG4caNG7Vs2TKVlJRo2LBhatSokWM/AAA1hdOu0Z8VFhZq+fLlKigo0NChQ9W0aVOzanOJrlEAQGW5HDUaGBioAQMG6OTJkxo+fLgZNQEAYJrfDMLc3Fz93//9n4YPH65bb71VX331lVl1AQBgCqdBOGXKFMXGxqpVq1b66quvFBcXp6CgIK1Zs8bM+gAA8Cinzwj79u2r22+//aqZZDIzM7VixQpTinOFZ4QAgMpyOmp04sSJGjBgwFX7v/32W48WBACAmZwG4a9D8OfJtvv27evZigAAMJHTZ4Rz5szR9OnTJV1ZRql79+7q06ePVq5caVpxAAB4mtMgPHHihKZMmSJJmjt3rmbMmKFNmzYpNTXVtOIAAPA0p12jzZs3l5+fn+x2u7Kzs9W3b1/5+PioWTPvDA4BAMATnLYIf14kd9WqVerdu7d8fHwkXVmfEACAmsJpi/DOO+/UH//4R/n7+2vhwoWSpBdeeEF16tQxrTgAADzN5VyjVRnvEQIAKstp12hBQYFOnjwpSTIMQ6tWrdKnn36qkpIS04oDAMDTnAbh+PHj9dlnn0mSFi5cqA8//FAHDhzQtGnTzKoNAACPc/qMsF27dnryySclScuXL1dycrKaNWummTNnmlYcAACe5rRFGBQUJEk6cOCAAgMDHa9N/LwfAICawGmLMDc3Vxs3btTSpUs1ZMgQSVdess/OzjatOAAAPM1pi/Cpp57Stm3b1KVLF8XGxkq60kU6ePBg04oDAMDTeH2iknh9AgCqt99coR4AgJqOIAQAWJrLIDx48KAZdQAA4BVOR43+7MUXX1TXrl0d2z4+PoqIiNDgwYNVq5bLwwEAqNJctghvvvlmlZaW6pZbblFpaamOHTumY8eO6eWXXzajPgAAPMplk65NmzYaO3asY/vtt9/W2LFj9c4773i0MAAAzOCyRbh7924dPHhQJSUlysrK0s6dO1VWVqbjx4+bUR8AAB7lskUYFxenCRMmKCcnRy1bttQLL7ygjIwMRUdHX/fF5syZI0k6efKk+vXrpx49eighIUHh4eHKz8/XpEmT5Ofnd/1/CwAA3OQyCO+44w6tWrXKsX306FE1adJE7dq1u64LnThxQps2bdLKlSt1/PhxjRkzRqdOnVJkZKRGjhypxMREpaSkqH///tf/twAAwE0ugzAvL0/ffPONzp8/L0nasWOHkpOTr/tCAQEBunTpksrKymS323XbbbcpLS1NI0aMkCRFR0crNTWVIAQAmMrlM8IJEyaorKxMERERioiIUEhIiFsXCg4O1tChQ/X8889r5syZuu+++2S32xUYGCjpSlDm5+c7PT4pKUk2m63cDwAAleWyRdihQwc99thjju1u3bq5daG9e/dq27ZtWrBggQoKCjR48GB16dJFhYWFkqSioiKFhYU5PT4+Pl7x8fHl9hGGAIDKchmEpaWlSkxMVPPmzeXj46OUlBS99dZb132hU6dOOYIuKChIvr6+6tq1q9LT03XXXXdp3759iomJuf6/AQAAleCya3Tnzp2qVauW8vLylJubq4KCArcu1KNHD0lXujhffvllPfvssxoyZIhycnI0b948FRcXq0+fPm6dGwAAd7lchmn//v3lXpXYu3evOnbs6PHCKoJlmAAAleW0a3Tnzp2688479dNPP+mnn35y7He3axQAgKrIaddoSkqKJOnzzz9Xbm6u48fdrlEAAKoil12jBw8eVKtWrSRJ586dU3Z2tjp16mRKca7QNQoAqCyXg2XWrVvn+O/z58/ro48+8mhBAACYyekzwoyMDMfPz1OslZWVOd77AwCgJnAahGfPnnU8E8zNzZUk+fr6atSoUWbVBgCAx7l8RpiXl6eIiAhJkmEY8vHxMaWwiuAZIQCgslw+I3zzzTcdXaMrVqzQ8uXLPV4UAABmcRmEkZGRuv/++yVJQ4cO1cmTJz1eFAAAZnE51+iFCxfKbTNYxroq0w1MVy6AqsplEDZq1EjDhw9XkyZNlJeXp/vuu8+MugAAMIXLIBwxYoRiYmKUlZWl2267zfFyPQAANYHTZ4QlJSWSpKNHjyogIECdOnVS3bp19eabb5pWHAAAnua0RTh27FgtWrRII0eOVNOmTSVdeX3i2LFjeuaZZ0wrEAAAT3IahIsWLZIk/etf/1KvXr0c+7///nvPVwUAgElcvj7xyxCUpOPHj3usGAAAzOa0RdilSxeFhISUm03GMAwVFhbqwQcfNK1AAAA8yWkQjho1SnFxcapVq/xH1q5d6/GiAAAwi9OuUV9fX9WqVUvJycnl9gcGBnq6JgAATOO0RXjgwAElJSVp+/btOn/+vGP/jh07rnpuCABAdeW0Rfj8888rMjJSwcHBioiIcPyEhISYWR8AAB7ltEUYHh6ue+65RzExMWrYsKFjf8uWLU0pDAAAM7icYu3SpUtKTk52dI/u2LHjqueGAABUVy7fI5wwYYLKysroGgUA1EguW4QdOnTQY4895tju1q2bRwsCAMBMLoOwtLRUiYmJat68uXx8fJSSkqK33nrLjNoAAPA4l12jO3fuVK1atZSXl6fc3FwVFBSYURcAAKZw2SKcPn26oqOjHdtHjrDSOACg5nDZIqxfv76efvppDRw4UM8888xVU64BAFCduUy1BQsW6KGHHlLTpk11+PBhJSUl6dVXX3XrYrt27VJqaqouXLigEydO6KWXXlJCQoLCw8OVn5+vSZMmyc/Pz61zAwDgDpdB2Lx5c3Xv3l2SFBkZqfT0dLcudOnSJX388ceaNWuWJCkrK0srV65UZGSkRo4cqcTERKWkpKh///5unR8AAHe4DMLs7Gxt2rRJERERys3N1aFDh9y60J49e1RWVqbFixfLbrfr/vvvV1pamkaMGCFJio6OVmpqqtMgTEpK0ty5c926Nqq3vIhmbh8bkcczbQC/zWUQPvXUU0pISFBWVpbatm2riRMnunWhEydOKCMjQwkJCSooKNCoUaMUHh7uWM0iICBA+fn5To+Pj49XfHx8uX02m82tWgAA+JnLIGzYsKESEhLk7++vwsJCt5dhCgwMVNu2beXn56f69evLMAz5+vqqsLBQklRUVKSwsDC3zo2Kq0zrCgBqIpejRseNG6cNGzZIkr777ju9++67bl2oQ4cOOnXqlCTp8uXLunjxovr16+d45rhv3z7FxMS4dW4AANzlMgg7duyoQYMGSZIGDRqk0tJSty7UoEEDDR48WHPmzNErr7yiSZMmaciQIcrJydG8efNUXFysPn36uHVuAADc5bJr9OdWnCQZhqGTJ0+6fbFhw4ZdtW/q1Klunw8AgMpyGYR33nmn+vXrp/r16+vMmTNuD5YBAKAqchmEAwYM0N13361Dhw6pefPmqlevnhl1AQBgigrNl1avXj117NjR07UAAGA6l4NlAACoyZwG4ZQpU1RaWqqzZ8+aWQ8AAKZyGoQtWrSQn5+flixZUm7/J5984vGiAAAwi9NnhDk5OfrLX/6io0ePatu2bZKuvD5x7NgxPfLII6YVCACAJzkNwtdee00nTpzQkiVL9Oc//1nSlSBcvny5acWh5mBqNwBVlY9hGIarD+Xn5+vw4cNq3ry5QkNDzairQmw2mzIzM71ag7d+wbu7qoLVAonVJwC44vL1ibVr12rOnDlq0KCBzpw5o3HjxjmmXAMAoLpzGYR79uzR119/LR8fH5WVlenVV18lCAEANYbL9whDQkLk4+Nz5cO+vgoJCfF4UQAAmMVli/DixYt69tlnHSvUR0ZGmlAWAADmcBmE48eP13fffaesrCzFxMSoR48eZtQFAIApKjTXaK9evdSrVy9P1wIAgOmYaxQAYGkEIQDA0lwGYb9+/ZSRkWFGLQAAmM5lEPbp00dt27Z1bB84cMCjBQEAYCaXg2XOnTunN954Qy1atJCPj49SUlL01ltvmVEbAAAe5zII09PTFRERoby8PElSQUGBx4sCAMAsLoNw+vTpio6OdmwfOcIkxgCAmsNlEIaFhWnChAkKDAxUjx491LhxYzVrZq0VDAAANZfLwTJvv/22hg8froYNG6pXr15asWKFGXUBAGAKly3CFi1aqGvXrtqzZ4/8/f3VqFEjM+oCbgh3119kHUPAOly2CLOysvTtt9/qzJkz2rp1q2PQDAAANYHLFuEzzzyjmTNnKisrSydPntSkSZPMqAsAAFO4DMLGjRtr9uzZOn36tBo0aCBfX2ZlAwDUHC5TbcOGDerZs6eGDh2qnj17asOGDZW64L59+9S+fXsVFhbKMAzNnj1bCxcu1LRp01RcXFypcwMAcL1cBuHq1au1fv16bdq0SevWravUqNGSkhKtXr1aDRs2lCRt3rxZly9fVlxcnNq3b69ly5a5fW4AANzhMgg7deqkoKAgSVJwcLDatWvn9sUWLVqkv/71r/Lx8ZEkpaWlKSoqSpIUHR2t1NRUt88NAIA7nD4jnDt3riRp//79OnjwoGOatdOnT7t1oX379qlevXrlXsa32+0KCAiQJAUEBMhutzs9PikpyVETAAA3itMgTE9PV79+/RQREeHYFxERoY0bN7p1oc2bN8vHx0fvvPOOzp07p8WLF6uwsFBFRUWSpKKiIoWGhjo9Pj4+XvHx8eX22Ww2t2oBAOBnToPwX//6l+NZ3i9FRka6daGxY8c6/nvp0qV69NFHtWvXLm3evFn33nuv9u/fr27durl1bgAA3OU0CH8Owby8PH3zzTc6f/68JGnHjh1KTk52+4LJyck6d+6c3nvvPcXFxSk1NVXz5s3TiRMnNGXKFLfPCwCAO3wMwzB+6wMPP/ywBg0apODgYEnSxo0bq8x6hDabTZmZmV6twd0pvCrL3SnAvFVvdcMUa4B1uHyhvl27doqNjXVsd+zY0aMFAQBgJpdB2KVLF02ZMsUxaKayXaMAAFQlLoMwOTlZf/rTnxxdowcOHPB4UXCNLk4AuDFcBuEdd9xB1ygAoMZyGYSnTp3Sc889p6ZNm0qiaxQAULO4nGLt6NGjiomJUUREhCIiIhQSEmJGXQAAmMJli3DWrFnlZpfp1auXRwsCAMBMLoNw+/bt2r59u2M7JSWlyrxHCABAZbkMws8//1wxMTGSpOPHj6tWLZeHAABQbbhMtRkzZpTrGv3oo488WhAAAGZyGYRHjx7V0aNHJUmFhYX6/vvvNWLECI8XBnhTZd7TZHo2oHpxGYSvvPKKY/HcwMDAcqtIAABQ3bkMwunTpys6OtqMWgAAMJ3T9wgXL14sSYQgAKBGc9oi/OKLL66aV3TLli0yDEObNm3yeGEAAJjBaRBOnDhRd999tyTp4sWLmj59ulq0aKHZs2ebVhwAAJ7mtGv05xDMzs7W8OHD1bBhQy1evFgNGjQwrTgAADztNwfLrFixQnPnztVLL72k7t27m1UTAJO5+7oIr4qgJnDaIpw8ebI+++wzffLJJ+VCkJUnAAA1idMWYVZWlmw2mxITE8vtz8zM1KhRozxdFwAApnAahP/4xz/0u9/97qr9W7Zs8WhBAACYyWkQXisEf2s/gMpjajfAfC4X5gUAoCYjCAEAlkYQAgAsjSAEAFgaQQgAsDSCEABgaS7XIwRwfSrzCgQA85kWhJmZmUpOTlarVq2Unp6u8ePH6+abb1ZCQoLCw8OVn5+vSZMmyc/Pz6ySAACQj2EYhhkXSk1NVXBwsNq3b6/169crNTVVUVFRunTpkkaOHKnExES1b99e/fv3r/A5bTabMjMzPVi1a/zrH1bGS/yoCUx7RtitWze1b99eklRWVqa6desqLS1NUVFRkqTo6GilpaWZVQ4AAJK8MFjGMAytX79eo0aNkt1uV2BgoCQpICBA+fn5To9LSkqSzWYr9wMAQGWZPlhm3rx5io2NVaNGjRQaGqrCwkJJUlFRkcLCwpweFx8fr/j4+HL7CEMAQGWZ2iJcsmSJoqOj1blzZ23cuFExMTFKT0+XJO3bt08xMTFmlgMAgHktwq1bt2r+/Plq06aN/v3vfys4OFiJiYmaOXOm5s2bp+LiYvXp08escgAAkGTiqFFPYNQo4F2MGkVNwMwyAABLIwgBAJZGEAIALI0gBABYGkEIALA0ghAAYGkEIQDA0ghCAIClEYQAAEsjCAEAlkYQAgAszfRlmABAqtw8vcxxihuJFiEAwNJoEQJwG6uvoCagRQgAsDRahACqHZ4v4kaiRQgAsDSCEABgaQQhAMDSCEIAgKUxWAYAKsjdQToM0KnaaBECACyNIAQAWBpBCACwNIIQAGBpBCEAwNIYNQrAUrwxUbjVpoSrbqNraRECACyNFiEA4CpWWmKrSgThqlWrlJOTozNnzujxxx9XZGSkt0sCAFiE14Pw3Llz+vTTT/Xxxx/r8OHDeu2117RgwQJvlwUA1Z6VWnWV4fUg3L17t1q2bClJuvXWW/Xjjz9e83NJSUmaO3fuVfttNptH63MpKMC71wdQs1Xmd1x1+/10A3+fZ2ZmVvizXg9Cu92ugIDy/7MuXbqkm266qdy++Ph4xcfHe6QGm812XTfNqrhPFcN9qhjuU8VwnyqmMvfJ66NGQ0NDVVRUVG7fr0MQAABP8XoQdurUSdnZ2ZKkw4cPq127dl6uCABgJV7vGg0JCdFDDz2kWbNmyW63a/Lkyd4uCQBgIX7Tpk2b5u0i2rZtq9///vfq27evQkNDvVJDTEyMV65b3XCfKob7VDHcp4rhPlWMu/fJxzAM4wbXAgBAteH1Z4QAAHgTQQgAsDSCEABgaQQhAMDSCEIAgKV5/T1Cb2Pli4rp06ePIiIiJEmjR49Wz549vVxR1bF3715NnTpV48aNU+/evVVSUqKEhASFh4crPz9fkyZNkp+fn7fL9Lpf36cVK1Zo6dKlql27tiTpgw8+8HKF3peZmank5GS1atVK6enpGj9+vG6++Wa+T79yrfuUlpbm9vfJ0kHIyhcV98ADD3hsrtfqLi8vT1FRUY7tlStXKjIyUiNHjlRiYqJSUlLUv39/L1ZYNfz6PknSG2+8oaZNm3qpoqonPz9fI0eOVPv27bV+/Xq9++67ioqK4vv0K9e6Tx07dnT7+2TprtGKrnwBaceOHXrvvfc0f/58nT9/3tvlVCkDBw4st52Wlub4hR8dHa20tDRvlFXl/Po+SdLHH3+shQsX6osvvvBCRVVPt27d1L59e0lSWVmZ6taty/fpGq51nyT3v0+WDkJnK1/gahMnTtQTTzyhdu3aacaMGd4up0qz2+0KDAyUJAUEBCg/P9/LFVVNXbp00ZgxYxQXF6cNGzZox44d3i6pyjAMQ+vXr9eoUaP4Pv2GX96nynyfLB2ErHxRcdHR0ZKuTJK+a9cuL1dTtYWGhqqwsFCSVFRUpLCwMC9XVDU1a9bMMaVip06dtHv3bi9XVHXMmzdPsbGxatSoEd+n3/DL+1SZ75Olg5CVLypm69at2rx5syTp6NGjjkEzuLaYmBilp6dLkvbt28c8kU7MmTNHpaWlkvhe/dKSJUsUHR2tzp07a+PGjXyfnPj1farM98nyc42uWrVKWVlZstvtGj16NKNGryEjI0Nz585VdHS0cnJyNHr0aLVu3drbZVUZa9as0eLFi9W6dWvFxsbKZrNp5syZql+/vgoKCjR58mTLj/KTrr5PP/zwg7Kzs9WwYUOdPn1aU6dOla+vpf9trq1bt+rZZ59VmzZtJEnBwcFKTEzk+/Qr17pPXbt2dfv7ZPkgBABYm7X/+QUAsDyCEABgaQQhAMDSCEIAgKURhAAASyMIgRtk27ZteuCBB/T3v/9d8+fP1/z58/XEE08oNze33OfWrl2rzp07O7ZHjRqlCxcuSJKSk5NdXue3jr+WipwTsDJenwBuoOeee0533XWXhg0bJkn64Ycf1KZNG4WEhJT7XJ8+fZSSknLV8c72u/u56/0sYEWWXn0C8KQVK1YoKipKISEhKisr09SpUxUUFKTw8HCVlJRIkr777ju9+OKLWrJkiTIyMnT27FklJSUpKipK/fr1c5yrIsfXrVtXM2bMUOvWrXXw4EH97W9/065duxzn7Nq1qxo3bqzZs2erQ4cOSk9P14QJExQaGqqnn35afn5+atu2rbZt26bRo0erb9++Onv2rGbNmqVbb71Vhw8fVpcuXXTvvffqgw8+0LFjx+Tv76/g4GA9/vjjXrnHwA1hALhhJk+ebIwaNcqYPn268eijjxo//vijYRiG8c033xjx8fGGYRhGYWGh0aFDB8cxI0eONI4cOWIYhmH07t37muetyPFff/21ERcXZxQXFxtHjhwxTpw4cdU58/LyjIyMDMMwDOPrr782ZsyYYRiGYaSmphrDhw83DMMw9u7da8TFxRmGYRgJCQnGhx9+aBiGYRQXFxtffvmlkZWVZdx///3lrn/w4EG37hdQFdAiBG6wQYMGadiwYTp58qRj1YDs7Gw1a9ZM0pUVBOrXr39d56zI8X/4wx+Uk5Oj2NhYNWvWTM8///xVn7npppu0evVq1atXT7m5ueVWW2nevLkkKSwszDHJ84EDB9S1a1dJUp06dTRo0CCtW7dOpaWleueddyRJjRs31pkzZxxLmgHVDYNlAA9p2LChcnJy9L///U8tW7bU4cOHJV1ZQeDMmTPXPMbX11eGYSgjI6Pc/oocn5WVpfvuu0+ffvqpbrnlFq1evVqS5OPjI+nKnLHvvvuu6tevr7i4OA0YMKDc8T9/7pdsNptjsM/58+e1du1a3XbbbQoKCtKYMWM0ZswYDRkyxBGiQHXkN23atGneLgKoCbZt26ZVq1bp9OnTOnTokLZv3661a9cqOjpaXbp00ffff68tW7bowIEDysjIUFBQkPLz87VmzRpdvnxZd999tw4dOqQtW7bo0KFD6tatm+PckZGRLo9v0qSJFi1apOzsbB05ckSPPPKI6tWrpx9++EH79++X3W5Xt27dtHTpUuXm5mrPnj1KT09XdHS0PvvsM2VkZOj222/XunXrHIvBDhgwQCtXrlRWVpZSUlLUu3dvtWrVShcuXNC6deu0Y8cOHTlyRH379r1mkALVAaNGAQCWRtcoAMDSCEIAgKURhAAASyMIAQCWRhACACyNIAQAWBpBCACwtP8HA8UTlXYNeAsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x311.496 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "width = 7\n",
    "height = width / 1.618    # golden ratio\n",
    "fig, ax = plt.subplots(figsize=(width, height))\n",
    "\n",
    "# ax.hist(psms_different['edit_dist_norm'], bins=np.arange(0, 1.05, 0.05))\n",
    "ax.hist(psms_different['edit_dist'], bins=np.arange(0, 25, 1))\n",
    "\n",
    "ax.set_xlabel('Edit distance')\n",
    "ax.set_ylabel('Number of conflicting SSMs')\n",
    "\n",
    "sns.despine()\n",
    "\n",
    "plt.savefig('iprg2012_consensus_distance.pdf', dpi=300, bbox_inches='tight')\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# conflicting PSMs with high sequence similarity: 275\n",
      "# conflicting PSMs with low sequence similarity: 467\n",
      "    (sequence similarity threshold = 3 amino acids)\n"
     ]
    }
   ],
   "source": [
    "threshold_different = 3\n",
    "num_high_sim = len(psms_different[psms_different['edit_dist'] <= threshold_different])\n",
    "num_low_sim = len(psms_different[psms_different['edit_dist'] > threshold_different])\n",
    "print(f'# conflicting PSMs with high sequence similarity: {num_high_sim}\\n'\n",
    "      f'# conflicting PSMs with low sequence similarity: {num_low_sim}\\n'\n",
    "      f'    (sequence similarity threshold = {threshold_different} amino acids)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEDCAYAAAAFhGKSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XdgHPWd9/H39l31XqxiSS5yr9jYgI0B06upCRAuhSSX5HJPyCXkyD2XSy7JJUf6c+RyIbkkEBIgoYRiMKYY2xjce5Mt27KaZXVpV9t35vlj5bWFbCODvcLm8/prtdN+O2vPZ38zv/mOxTRNExERkSSxDncDRETko0XBIyIiSaXgERGRpFLwiIhIUil4REQkqRQ8IiKSVAoeERFJKgWPiIgklYJHRESSSsEjIiJJpeAREZGksp/qAtXV1WeiHSIicparqakZ0nynHDwAK6tGvZ/FRD7Upr+9lKYHpw13M0TOStU/9Q55Xp1qExGRpFLwiIhIUil4REQkqRQ8IiKSVAoeERFJKgWPiIgklYJHRESSSsEjIiJJpeAREZGkUvCIiEhSKXhERCSpFDwiIpJUCh4REUkqBY+IiCSVgkdERJJKwSMiIkml4BERkaRS8IiISFIpeEREJKkUPCIiklQKHhERSSoFj4iIJJWCR0REkkrBIyIiSaXgERGRpFLwiIhIUil4REQkqRQ8IiKSVAoeERFJKgWPiIgklYJHRESSSsEjIiJJpeAREZGkUvCIiEhSKXhERCSpFDwiIpJUCh4REUkqBY+IiCSVgkdERJJKwSMiIkml4BERkaRS8IiISFIpeEREJKkUPCIiklQKHhERSSoFj4iIJJWCR0REkkrBIyIiSaXgERGRpFLwiIhIUil4REQkqRQ8IiKSVAoeERFJKgWPiIgklYJHRESSSsEjIiJJpeAREZGkUvCIiEhSKXhERCSpFDwiIpJUCh4REUkq+3A3QEROj396KcBXLnRRlnlmfk8eWb8F+M/lIep7DK4cY6e1zyTTZeG+i1w8ujHM998M8eW5TkwTth2O8Q9zXUwpsgGwvjHKa/uiZLkt+CNwoNPg7893MrHQxveXBUlzWrBZob7b4DsL3XgcFvwRkwdXhChOt9DQbXLfRU5yU+KfceneCP+5IsR3F7q5YGT8cLb1UIzHNoepzrexpz3G5aPtLBztOCP7RN4fBY/IOeIn13iStv5LRtl5+2CUL85xAXD9I31cXW3nnhlOvv9miL8/34ndauHVvREeXhvmoRs87GmP8as1YX5zswerxQLASzURGnoMJhbayPZYEuv79utBnt4e4e7pTv6wIcyYXCt3TXOyfH+Un70V5ntXuGnvM0h1WhiRPjBoOwImn5rpZHyBjQ6/wTV/8Ct4PmQUPCIfQn/dFuZXa8K8cW8ar9VG+OdXgqz/UjrP7Yzw4IoQn5zpoK7LpMVr8PAiD1taDL7zepBvLnBxfpmdP20Os7Epxug8K8v3R6nItnLLJAffWxbimwtclGRYeeCVIIsmOrh5ooOathh/2BhmVI6Nui6Dr81zkeWxJNqzsTk2YP3HCsdM/FGTdJfl3R+DroBJTkr8/ce3RLim2pEIHYBrqo8GwpHQATBMSHXG51tZF+OBBfFpU4tt/OtrQQDyUq3kpVr5nzXhAdu8pMr+rvUcfx//fFUImwWiBgSjJg8scPNOfZSX90QpzbBwyGvywAIXL9dEj7vPIzH4zhtByrOsNPea3Djeznmldv64Kcwhr4HTZiHdZeEz5zn52aoQz+2McPc0J+ubolTn2bjvIhf7OmP8ek2YMXlWdhyO9/J+tyH+3f3x9hSe2RHhoXdCvHFv2nHnzXQP3udnA13jEfkQum3y0aPlwtEOMvoP6jdOcFCVY2VCgY3vX+HGboVdbQYzRtgYnx//79zpN/jNujAPXu3mC+e7sFotLJroYFapPTFPaaaV2aW2xDb+76tB7pnu5N5ZTmaU2PjdhoEH82PXf0Rth8HDa0P87K0Q//cSNyUZR6f/dl2Y+18OsLohxj9fHA+Npl6D3JSjB8oleyJ8b1mQh9eGBqz3sNegqcfg2nHxAOkMGInwSHNBZ8Ac8n58ZGOEBy52DXp/+YEoB7oMvnyBi/suclGdb8M0Tb6xJMjX57n43GwXaU4LT26NnHCf7+8y2NNucNdUJ1+9yEVuqpXajhjP7Ihw/3w3X7nQxZv7o+zvNLjvQhedfpOPTXXwyxs8PLMjAsDKAzE8Dgufmunk78934rDBbZOOhvHNE4++Pt68ZysFj8hZqDwr/l8322OhLzzwQNzQY5KfasFmjR/kSzPe+1fx3naDZfujPLw2RG1HDOsQfkiPzrXyudkuvnGxe0AvA+DeWU6+d4Wb2g6DQ14DgJIMK219RmKeq8Y6mFBgY0/70fc6/QY/WhniR9e4cdrijcjxWOnrz0FfCHI8Q/uV/+ctYcozLVw+ZvBptj3tMUqO2S83T3TQGTCJxEj03EozLew9pm3v3ucTCmzcOsnBF54L8MArAewW2NthEDPg4bUhHl4boijdQqc/vo7cFEv/NSxLIjRumxz/UXH7436e2Bo56QH5VOb9sDub2y7ykRCOmnT6B4bLyQ69ZZkWDvtMYkZ8mcbeo8umOi2Jg3iL7+j71flWrql28LnZLu6d5eS8kg/+c9pps3DvLCe/Wh3f4J1THby4O0rUOLrd2DGv2/sMfrg83nvK8Vh5fV8UgHkVNra1xADYfCjG/Ir3vkLwp81h3HYLt09x8vbBKMHIwP1XnWejsefoe3/dFibHY8Flh95g/P2GHpOxeUcPke/e5w09BjNLbPzpjhQurrTzyMYwY3OtpLksfG52vNd080QHI7Pj67Ac50vb2hLjC3OcPH1nCr1BkxV10f7vKN6GFq9x0nnPVrZvf/vb3z6VBR566CE+nZ1zhpojMnx+3bCPf7q8aLibkdDeZ7Jsf5S2PpNNhwzSXeANmfxtZwSrFTJcFp7cFqHDDzkpFp7ZEX99xRg7bruF32+M0NBt0NRrMqvURmmmlUyXhUc2hekKmBzsNtjfaTC71MYF5XZ+vyHC7rYYr9XGuGSUfcD1g83NMZ7uX39ljpXHt0Q40GWQ6bYwJu9oSD26MczKuhgOK8wqtTM618ov3g7T6Te5ptpBWaaF32+IsK0lxpsHouxoNbhrmoMRGVbufNJPh9/k9X1Rnt0RIWKYXFRhZ2Khjed3RdnVGmNTs8FX5zlJccTb9vDaEO/Ux+iLmOSkWClOt/JabYT/WBairc/k2R0RXtgd4dZJTlz2o59nZLaV3e0GqxtirG6IUZ5pZVSujfEFVn6/IcKuthjdQZO/P9/JusbYcff5tGIbv1wdpr7bYG+HyS0THYwvsBGMmrxcE2F9U4yGHpPLRtl5dmeU1/dFGZ1r40CXweKaCAVpVqIGPLE1Qm2HQXfQ5GNTnGR7LGxsjrHjsIEvDK/VRinLshKMDp7X7fjwXON56J0wX/7yl4c0r8U0zaGfMAWqq6tZWTXqfTVM5MNs+ttLaXpw2nA347T75yUBFk10DBoUIHI6Vf/US01NzZDm1ak2kXPYtpYYNe0Gz+0cfLpJZLjoJ5DIOWxykY1n7049res0TRMTBgyLPtcYpokFsJzDn3E4KXhEZMi8IZOfrwrxtXku1jZEWVobpTLbyq7WGHdNdzJjRHxY8oMrQqQ6LXQHTM4rtXHV2PjIssYeg/94M0SKA358zA2p33sjSLrLgkn8ov0PrnDjtA886K9tiPLSnigFqRZafCaRmMkProyv46F3QkQMcNpg+2GDG8bbuXqsg/9dH+bBFSH+4wo3txwzTPmtuiifeSbA1+Y5+eyswcOtQ1H4yVshvnKBi7Tj3J8kH4xOtYnIkH1vWZBbJznwOCwc9hncd2H83p+/m+Hke8viN3a+VhvlsM/kH+a6eGCBi5+sDNHdf+/N1pYY8ysHj5gbmW3l/1zo4isXughEYNn+wSO2frg8xNcucvHFOS7+faEbR/9w65q2GLvb4vfKfGmOi69c4KSjL769z5znZGKBlYfXhQaMoHvrYHz9n5p5/LtLPQ4Lt01y8N3+zySnl4JHRIbEFzZZ1xhjfEE8OG6f4iQvNX4IMUxI7R9hdbDbpCg9/tpmtZDhtrC1fzj0NdWOxP05x/rE9HgAxAyTVp9BRfbgQ1NuioVfrQlz2BcfYvzvC90AZLktbDkU49W9EcJRk+p8G3dPPxooY/OsVOXYeKkmHjbrG6MDhoubpsm/vRbkoXdCfG9ZkCV74jd3Vufb2NgUwxvStbHTTcEjIkPS2GOQdYISLX/eEubr8+KnrKaNsLKtJR4O3pBJXZeBL/zeB+/1jVH+4YUg8yvtVOcP7hX959VuwjGT2//s5+Y/9SUCojDdyk+vdfP0jgjzHvbxlRcDNPYYA5b9wmwn/7M2jGmavLYvymWjjl5l6AnCigNRbp/s4BvzXQMqMGS6LdR3D1yXfHAKHhEZknAU7Mc5YvxiVYhrqh1MKY6HxXkldu6Y4uCXq0M8vT3C5CLboEKex3NeqZ1f3eihM2Dy6MbwoOk5Hiv/comb5Z9L44GLXfzba0Ea+kNhVqmd/7kphWWfTaMiO16H7lhTim0UpFr46VthJhbYBgwayPJYuP9iF/cvCfLZZwMDbvS02yAUG8rekVOh4BGRISlMt9DzrtNOP18VYkaJjQVV9kSlgUAkXk7mS3Nc3D09fkF/ctHJDzVHapcBlGVaqe8Z3Mv411eDHLntcFapnYJUKxEj3hP779Xxem8pDgsLR9sJxwb3sL5wvpOltRGuqR44pqo3GC8x9IdbU/j0TCcPvRM6ZhoUp2twwemmUW0iMiSFaVYyXPGRalkeC49sDPPU9ggbmmI8vDZMZyB+l35f2ORfXw1ywUgb/jD8ywJXom7c8gNRlu2LUt9j8PiWMB+fGr8Ws+pglENeA9OEmjaDb14yeKTZiAwL33otREFqvCTQtePsVOVY6Qma7G4z+Olb8WrT+zoNHrg4fv3n0Y1hatoN1jVGmV1m55VPpQHw2KZ4j+qPmyLcNMHBoxsjbGyK0ew1uWNKvE29QZM0JxQPobcmp0aVC0T6nauVC06nbS0xnt8V4V8ucQ93U864H7wZ5LpxDiYXncVloJPoVCoXqMcjIkM2uciWqM585Hk556K+sMknpjspPUNPc/2oU/CIyCn5KByMU52WczpYh9u5/y9IREQ+VBQ8IiKSVAoeERFJKgWPiIgklYJHRESSSsEjIiJJpeAREZGkUvCIiEhSKXhERCSpFDwiIpJUCh4REUkqBY+IiCSVgkdERJJKwSMiIkml4BERkaRS8IiISFIpeEREJKkUPCIiklQKHhERSSoFj4iIJJWCR0REkkrBIyIiSaXgERGRpFLwiIhIUil4REQkqRQ8p5FhmnypuWnAe2/4fFxdd+B9r/MvPd2J192xGA+0HDrldbzq836gNoiInE4KntPIarHwyxElA967NC2NNOv7381/7elJvM6y2fhBUfEpr+PytPQP1AYRkdPJPtwNOJe86vPy0/Z2Fo+s4MH2NlKtVnJsNiKmmZjnqZ4eWqNRnBYLqVYrH8/K4jedHSzx+rglM4OtwSBVTiefy8nlJW8vPsPgd52dTPd4qAuH+XNPN38tH0nUNPmvjnbybHa6YjHy7XY+npXFD9taKbLb6YjFmOhyc1V6+jDuERGRwfQz+DQ60rN4y++nzzD4cm4eN2dk4jUMAA6Ew7zk7eWLubncm5PDO/4+6sNhPpuTS7cR46aMTL5fWMTLXi8A16RnkGa18umcHKZ7PCzKzExs6wVvLy6LlU9kZ/OPeXlk2eJf5byUVD6ZncN9uXk80t2V/J0gIvIe1OM5AxoiYUbYHQB4rFayrDYA6sJhDOCx/kDIt8d7K+VAts1GSv/pMLvF8p7b2B8OU+VwJv6+Oj0DgMPRKL/p7MBjtdIdi53GTyUicnooeM6AMoeTpcF4ryVgGHQb8QCodDpJtVq5OysbgPUBP6WOeECdKGqOvF8bCjHa5Uq8P8rppCkSAcA0TV7yeql2uXi2t4c/lpUD8PQx14dERD4sFDyn0es+Hz7DoDsWI8Vq5WftbeTY7KRarTzf28sNGRksTE3jlx3tuCxWYpjMdHsS13JW9fUB4DMMlni9XJWezliXi191dJBnt7EzFMJnGCzz+bguPYP/6mjnD12d9BkG81JTKXM4KLI7+HFbG8UOO0HTZInXi8NiwWcYvNDby/UZGcO8l0Tko85imsdc+R6C6upqVlaNOlPtERk2099eStOD04a7GSJnpeqfeqmpqRnSvBpcICIiSaXgERGRpFLwDJFpmhindlbyrBT9CHxGERleGlwwBD4jxm86O/lCTi7u/qHOW4MB/rG5mcfKyhMj0w5FIvy/jnY8VivfKihMLL8+4Gdxr5dKp5M94RDfyMsn3Wbjd52dbAoGEvN9PS+fcqdzwLb3h0M81dNDkd1BjxHjYDjMj4tHAPBsTw8t0SgpViv7wiGmuj3ckpnJC729/HdnB7dmZPKZnJzEug6Gw9zT2MCtmZl8PicX53GGbf+6s4NbM7MotOufhoicGTq6DMEv2tu5PTMLd/99Nt5YjFd9PvLfdXDeFQoxJyWFLcHggPf/p6OTr+XnMc7l5hft7SzxebktMwuA/3pXiZ13+11nF5/Mzk4Mpf5RWxsQD8PFXi+/LS0FoCcW48n+um7XZ2Swoq+PF729fCwri9T+di/xecm22ViUkXnc0AG4MyuL77a28tP+cBMROd10qu09+A2DzcEgY465h+a3XZ18Ojt70LyXpqXhOM4BPdtmoycWr17QZxhUHdOrebSriz91d/FUTw+x45zmyrbbeLKnJ3HPztfz8wFwWqy0x6K87O3Fbxhk2mx8Lic3sVyWzcaC1LTEvTzNkQh5NvuA9v2ms4PfdHbw684OnuqfL9tmJ2SaHAyHh76TREROgYLnPTRHImQcU2Dzhd5eLkhJJds29M7i53Ny+HN3N7/saKfXiDHKGQ+xS9LSuD0zk7uysmmPRfnLcW74/GJOLoV2O/90qJm7G+p5tjc+j9Ni4cdFxaz2B7i9/iBfP3SIvaHQgGU/npXFs709BAyDxd5ern1X3bYX++8V+nxOLtXHBGuuzUZdRMEjImeGguc9hE1zQAmb7cEge8MhHuvuwmcYPNfbw453nVo7VsQ0+UZLC98uLOBLuXnM8qTw3x0dQLySwZHTd9PdngHXe47wWK3cm5PDE+Uj+beCQh7v7mZTID7faJeL7xQW8uzICi5MTeGbh1sGLFtgtzMnJYXfdnWSbrUltnXEvxYU8FBHB19saiLYX08O4iV7whpkICJniILnPRTY7YkinwAPFBRwd1Y2d2dlk2a1cmNGJhPd7hMuHzZN+gyD9P56bcdWqz4SQABNkQil/fXdjvWz9jZC/dsf43JR6XQmQuGHba0AOCwWLklNO25Y3J2VzVKvl+syBlepDhkm/1lUzDfy83mo82hbvDGDglPo0YmInAodXd5Dnt1OmtVKbyxGhs2WeP+RrniP59neHm7LyKTI4WC1v4+3+/w0RSP8rbeHmzIySbVa+WxODj9ub2OE3cH+cJjP5MSvD1mJD1zIstmoj4T5x9y8Qdsf7XTx4/Y2iu0Oeg2DEruDWR4PEK8D98uOdtwWK/WRCPfnxa//LPF62RcOsaqvjwtTU3mhojLx/pE2fzEnl5e8XvaFw/QaMW7NiFe+jpkm7bHoScNUROSDUMmcIdgdCrLU6+Mf8wYHw7nmT91dVDicXJiaOtxNSTqVzBF5/06lZI56PEMwzuUm02rDbxiJRxecq+ampFDldL33jCIi75OCZ4iKHYOvv5yLFDoicqad2z/fRUTkQ0fBIyIiSaXgERGRpFLwiIhIUil4REQkqRQ8IiKSVAoeERFJKgWPiIgklYJHRESSSsEjIiJJpeAREZGkUvCIiEhSKXhERCSpFDwiIpJUCh4REUkqPY9HpF8oCiPu3zrczRA5K6UXVQ55XgWPSD+TDP7rM48NdzNEzkrfXHzfkOfVqTYREUkqBY+IiCSVgkdERJJKwSMiIkml4BERkaRS8IiISFIpeEREJKkUPCIiklQKHhERSSoFj4iIJJWCR0REkkrBIyIiSaXgERGRpFLwiIhIUil4REQkqRQ8IiKSVHoQnIgM2ZOb/sjl1deQk5J7xrZhmAY/euO7LJp8O2MLxgPQ0H2QV3a/iGnCZ+Z8AavFSmN3PUt2v8CIjFKuqL6WpTWLWXNwFV+46D6K0ouJGTFe2f0izb2NXDnuOsqyRg7Yzuq6t+gOdOG0u2jxNlOZM4q5FfOO2ybTNHl51/Osa1jNDRNvZnrprDP2+T8K1OMRkSG7Y/onzmjoANS272FkTiXrG9ck3ivLGklFzihiZpQ39r4CQGlWORX9YWG32ZlbMY/ijBKe2Pgo4VgYm9XG3Ip5VOSMGhQ6wUiADY1ruGr89Vw65gpumnQb3mDvCdtksViYWzEPl92l0DkN1OMR+YhaV7+aZbVLuf/Sb7GzZRtPbfkz37ryB2xqXM/Lu5/nosqLae9royfYw9/N+iwNXQd5fsdTXDthEVW5o1ld9xZ1XQcoTCuipnUneWn5zCw9nxd3PsO1ExaR7cnhqS1/ZkbpbGaWzaalt5lVB5aTn1ZIe18bV427jhRn6qB2Heio5ebJd/DjN7+PP9w3YJ6bJt3G45sepTJnNKPyxgxadnR+NTEjyvPbn+bWqR8/4We3Wx30BnvZ2LiWiUVTSXGmcsW4awEIRPy8vOt58lLzae9rZ27FRRRnlJxwXVEjyuIdfyPTk0V3oIvJxdOO2zY5Sj0ekY+oWeVzEq8nFE3G7fAAML30PPJTCyjOKOXmKR/DarFyqLeJkTmViQOwL+Rj+b7XuX3aXVwy5nIsFgszSmdTmTsqMU92Sg6VuaMT23hm25PMrZjP/FGXMjK7kpX73xzUpkDET6ozDafdxeTiaWxu2jBgusvu5mPT7+HZbU/iC3mP+7kur76GTn8HmxrXnfCz2212Pjn7c9S07uZHb3yXR9Y+THNPEwBv1r7GiIxS5o+6jHlVl/D0lidOuh/X16/GaXeyYPRCrh5/PX/d8icM0zjpMh91Ch4ROa7c1DwAUp2phKKhAdO6/B2kuzOwWuKHkKGcfjvsbWF36w6W175Gq68Fq8UyaJ6tzZvoCXbz2p4lGIbBhsa1g+YpzhjBgtGX89fNfwLMQdOtFisfm/4JXt+7lLa+1sT7r+15md+v/TWv7Xm5fz0lfHzGPfzzwm8zrnAij23430Q7s1NygHh4HvYdes/Ple2Jz++yuzFN84ShKHEKHhEhGovSF+4b8J6FwcFwRHZKLr3BnsQv+05/R2Ka0+YiFA0C0BvsTrxflF7MlBHTuXj0QuaPupSKnKpB6+0KdHLthJtYOPYqbpx8K1EjmuiJHOu8svNJd2ewvmHNoGkAGe5MFk2+nWe3Ppl4b+HYq/nU7M+zcOzVADyzNd6TsVvtTC6eRsyIxtuZUUynvzPxuYrSi0+4H1YdWE5RRjFdgfj8wUgQq8VKmiv9hMuIrvGIfKRNHTGDF3c8S2F6EW67m3X1q8lJyaU70MWmpvWML5xES28zmxrXY7faE69vnHQrF4+6jCc2PUpJZhk269FDybSSmbxZ+ypd/k5C0RC7W3cwtmA8N0+5gxX73iDLk01PsJv5VZcOaMuqA8tp7mlMXNfp8nfisrtYsvt5Lq++hrrOfYSjIa4cdx02q40bJt3Kr1b9HIgH5zt1K2nubaSxu57SrHJG5Y1hVvncE572CkVDvLzreRw2B22+NhZNvgOABaMWsnjncyzb+yqd/nZunvIxTNPknbqVhKIh3ti7NLGOTn87iybfwYs7n+WNvUvpCXRx69Q7Ez1BOT6LaZqD+6onUV1dzcqqUWeqPSLDZvyK1Tz8+ceGuxlnpSODCKqOuaYjHy3fXHwfNTU1Q5pXsSwiH0hjd32iJxSJhYe7OXIW0Kk2EflASrPK+Yd5XxvuZshZRD0eETklpml+JIYLx4zYcDfhnKUej4gMWTASYGnNS1w1/nqcNif7O2pZvPNZ5oycN+C+oN2Hd7C3vQaHzUFuSh6zyucCsKVpI7XtNWR6sunyd7Boyh3Y+wcmmKbJI+seJsWRyu3T7x607RbvId45sIIsTw7+SB9tvsN8cvbngROXv1lXv5qXdz3PBZXzEqPZAFp9h/nF8v/kgsr5XFl9HXbb4EPhK7tf5ILK+WR5sk/rPhT1eETkFLyw4xnOK5uD0+YkZsTwhnoH3dV/ZLTYtRNu4qpx17OhYS1d/cOTX9jxDFeNv56FY68iGA2yp3VXYrl36laSeZKD/Ot7ljCnYh6XjLmcayfclJj3ZOVvZpXPYWR2Bevr1xCMBBPr2tS4njRXOnNGXnjc0AGYP+rSAcOx5fRR8IjIkISiQQ507mNEZjxobFYbU0fMGDRfQ3cdOSm5iSHFJVll1LbvASDNlYY/7AcgGouQl1YAQFNPAxEjQnlWxQm3n+ZMZ9WBN+noawdg0eTbgYHlb0LR0IDyNwCprjQmFk9h9cGVQPzenAx3BjarLTHP0pqXWFrzEq/sXszbB1b0tzWdcCxMq+/wqe8sOSkFj4gMSae/gxTH4Npq79YX6sNpdyf+dtnd9IV9AFw7YREv7nyW57Y9Rbo7kyxPNuFoiNV1q5hXdclJ13v1+OvJdGfz+7W/5mfLf8jqg6uAk5e/OWJ+1aWsrltFOBpifcMaZpadP2D6hobVzCidxZXjrqUksyzxfro7g1avgud00zUeERmSqBEd0Es4kVRXKuHo0dNaoWiQbE8OvcEentv+V/7P/Ptx2Jy8tPM5Vu57g5LMMtwONyv3vUFTTyMd/nZW7lvGvFEDg8hpd3F59dVcXn01zT1N/GnD7yhIK6Qqd3Si/E3UiLKhYQ2Pbfhf7r/0W4llMz1ZjC0Yz6t7XibDnYnT5hyw7tum3c1LO/9GIBLg0jFXJt63WWxEjcj73WVyAurxiMiQZLqzCET87zlfWVYFnf6OxMi3xu56RueNxR/xY7PYcPQf9NNdGUSNGOMKJ3LthJu4ePRCqgsmUJhWNCh0AJ7f/jSRWDwERmQLRI5DAAAVMElEQVSWUJheRLS/zM2Jyt8c6+JRl7G5aQOzyuYMmhaJhbln1mdZNOUOXt71XOL9QCRApifrPT+znBr1eERkSDLcmbgdngGPKlhzcBUtvc0EI0Ey3JlUF4zHZXdx9fgbeGHHMzisDs4rOz9RdHNi0RSe3/40aa50Wn2HuW7CosT6D3TsY3frDjr97WxqXDfouTfFGSX8bdtfyE7JJRDxk5OSx+i8scCJy99sbFxHS28zuw5vZ3zhJP7l8u8CsKlxHcFIkNUHV3H1+BvY0LCWw95D+MN+5lbMB+IPpPOGek563UneH5XMEemnkjnvrbG7ns1NG7hu4qL3nvkst2Lf6+SnFTK+cNJwN+WscColc9TjEZEhK80qJ9WZRigawmV3DXdzzqixBRNOWpla3j8Fj4ickiOnzc51Cp0zR4MLREQkqRQ8IiKSVAoeERFJKgWPiIgklYJHRESSSsEjIiJJpeAREZGkUvCIiEhSKXhERCSpFDwiIpJUCh4REUkqBY+IiCSVgkdERJJKwSMiIkml4BERkaTS83hE+oWNKPc+/InhbsZZxTBiWK36/SpQUJI/5HkVPCL9DKeDGx98cbibcVb5679czuef+sxwN0M+BBZ/fsmQ59VPFRERSSoFj4iIJJWCR2SY9dTvpPaV/2X/64+y5Y/fYsdffgBAJOBl93O/OOmyrTtWEg32Jf7e+dSDxCLhM9rej5r23R1seXQb2/+8k5XfW8U7P1kDQNgXZv1/bzzpsg1vNxLxRxJ/r/7ZOmLh2Blt79lA13hEhln77tXkj59LRtl4gt2t7H/1dwA4POmMu/H/nHTZth1vkVE6Drs7FYAJt95/xtv7UdO87hAls0eQW52Dv83P1j9uB8CZ5uS8L8446bKN7zSROzYHR4oDgDn3zTrj7T0bKHhEPqD6VU9jsVgxjRhGNEJ21VQ69qzFlZFHyNtJxYI76ahZy8EVT1A88yqCXS2EvJ2MX/RVfC378B0+gGnEiIb8pBaMTKy3ZfPrNK9/iRn3/gQjFuXg8sdxpmUT7uvBlZ5DamEFod52Dm1cSmpBOQ5PBgfeeJQJt/0z3qa9x92exWqlZfPr9DbtJSVvBF37t+LJLqTysnuwOVzDuBfPnJrn9mKxghEzMcIxCqYUcGhDC55cD8HOIBNuH0fz+hZ2P72byoUV9B32E+wKMuvLM+k60E1PXQ9mzCDij5BZnpFYb92b9exfup9L/2MBRsxg5192485yEeoN4852kzUyg0BHgAOvHySjLB1Xuovtj+9gzldn01nbfdztWawW6t6sp6u2i7QRabRuayOtMIVJd07E5rQN4148vXSqTeQD6DqwlWBXC2UXLKL8oltJySuldslvKJ93OyWzr8Pm9HB46zLyJ1yAJ6eY1IJyRl3xaSxWG31t9aSPGENqfjn5Ey8id8x5A9ZdNO2yxOvW7Suw2p2UzL6Oykvuwu5JJ7NsPK6MPIpnXEHBxHlkV03FlZEHcMLtRfy9NK1bzJirP0fp+TdgsVrJn3jRORs6rdvb6Dvcx9jrxzDuprGkl6Sz5fdbGXfzWEZfVYXdY+PgigZK54wgtSiVjLIMptwzCYvNQm9DLzmjsskoS6d0bgnFM4oGrLtiQXnidcNbjdicNkZfPYqJd4zHle4ktzoXT66HystGUnZBKQWT8/HkegBOuL2QN8y+JfuZ9ukpjLlmFBarhdK5JedU6ICCR+QD8bc3Jg72ANlVUzFiUeyuFADcmfn425sS091ZhQA4PGnEwoH3vZ2CiRcNabl3by/Y04YzNQtL/703x67zXORt8iYO9gAFk/MxYiYOT/zUV0puCr5mb2J6akH8e3OmOYgGo0PeTu+7tlM6t2RIy717e/52P+5MFxarJd6+PM/JFj9rKXhEPoCUvFKCPW2JvztrN2K1OxIX/IM9raTklZ7W7ZimSev2FfEJlvh/4b62hiGtx52ZT9jXhWkYAIR62z9w2z7M0kvS8bf7E38f3tKK1W5NXPD3t/tJL0n/wNvJKEkn0L8d0zRpWNUIgMUSD5DeRu8Jlz1WSl4Kwe4gpmH2t2/oP07OJrrGI/IBZFdOwdu8l/pVTwOQWlDB6Ks+y8EVT+LKyCUWDlI4ZQE99bsI9XbQvuttskfNwN/eSNvOt7HaHInXDk867btXE+rtoPvgdkI97URDATr2rKVg0nzq3nycxtXPEQ0FyBk9E4Csisk0rX0Ri9VOzqhphHo7aN2+gsyyCcfd3qgrPs2IWdey96VfkVpQgcVqAyzDuAfPrIJJ+XTVdlHztz0AZI7MZNqnJrPr6Ro8OW4iwShj55XRXtNBoCNI05pmCqcW0tvoo3FNM1a7NfHameakae0hAh1B2na2E2gPEAlEObShhbKLStn5l93sXVxLJBClaFq8p5k/MZ/al/djtVspnJJPoCNIw1uN5I7PPe72pnxiEqOuqmLTb7eQWZ6B1WYBy7n3/VhM0zRPZYHq6mpWVo06U+0RGTaVb6/hto9A5QLvoX2kF8f/D+986kdUXf5J3JlDL3dyLFUuOP26DnSTXZkFwJqfr2Py3RNJyUsZ5la9t8WfX0JNTc2Q5lWPR+Qjpm3nKrrrtmMaUTLLx7/v0JEzo+mdZtp3tGPETHLH5Z4VoXOqFDwiHzFVl90z3E2Qk5h054ThbsIZp8EFIkk2lIoEJ5q/ZfPrbPztP52ppn3kNK8/xCtfeQ0YWFWgaW0z2x/fyebfb6W30cvGhzdT+9I+dv5l1xlrS8vmwwMGQpzLFDwiSTaUigQnmv/Ye3vkgxtxXjGOlPiJnzn3zUrcL9P4dhNVCyuY+neTiQajWB1WRl8zinG3VJ+xthzefJhAx7k5iu3ddKpN5AwLdDbTuOZFUvJK6TtcR0p+Ka3bljPj3p9Qv+pp2nauomjaQnoO7iCrcjLRgA9fywFyxsykaOqlAyoYHCsa8lP78sOkl4zF395A8YwrSSusYO9LvybQdYjM8gn01O+i9PzrEqPgBEzDZNtjO7C7bbgyXRhRk9ZtbYmqAr5DffS19nHgjYPkT8jj0IZDeBu97HlhL2OuHc2BNw4S7ApitVuxe+yMuqKSmr/toXF1M2UXltBZ203h5Hxyx+dy4LU60orS6GvtY9zNY+k+0MOOJ3ZRPLOIaCBKb2Mvs748E397gN4GL4000dvgpXJhxXDvpjNKPR6RM6zrwDZsDhcjZl5JyfnXUzzjysS08gtvIeL3UjTtMsZe/yUaVj1DyfnXM/b6L9Gy+Q3gxL0ci8VC8cyrKJl1DSNmXkXT2viIvLILbybS10P5hbcwftF9pBZWnfkPeRY5vLWVSCDChNvHM3JBORF/ZEBVgSOvKy8bScHkfErnlpBRls7Y68fga+mj8e0mxt9STfWNY2jd1oavxUf1TWMJe8NUXlbBeV+cTv7kfLY9up2KS0cy6spKskdlsX9pHQWT8skZnUVKnodJd04goyyD9l0dZJZlJCoknOuhA+rxiJxxhZMvpnHNC2x7/LukFVUxcv7HBkx3pmYkStY4UtITr43Ie512sdBzcDve5j3EQkGigaM3KboyC7BYbThSMk6y/EdTX0tfYqSY3WXHmeYc8rK+Qz5M06R2yX4APNluwt4wFIErw5koBppWaMfb7KN1axttO9qJ9EWw2o/+zk95nxUSzhUKHpEzzNeyn9I5N1B+0a3sfel/6K7belrW27p9BdGQn6qLbiXQ1YKvZV9i2jl4z+Fpk1qUSvfqZgCioShh39AfI5E2Ig27287oq+K9yPZdHaQUxCuDv/s+3PSSNEacV0RqYSphX5ie+t7EtON+PVYLJuBr8eHJTcHmOHdPSCl4RM6waNDHweWP48rIw2K1EexuTVQkiEVCREMBuvZvJhryD3rdsXcD0YA3MX806E+8zhw5iY696zm48i9gQqi3g56G3fQ2xKsktO1cRf6EC4f743/oFE4p4PDmVrY/vhNXhgu7287BFfWJqgI5Y7IJdASpe7OeyktH0rimmd5GH01rmik5fwQjZhWz66ndWB02TMMkd1xO/Lk7gSh1yw5ScUm8wviUeyaxb+kBPDkegl0Bqq6soqe+J1GlwJOXQufeLhxNPoqmFpI3LpeGlY1YbBam/t3kYd5LZ5YqF4j0+6hULjidVLlAjjiVygXnbl9OREQ+lBQ8IiKSVAoekffJNE1M0xjuZpwWphEb7iacdqZpJh4vcLYzYufGv7MjNLhA5H2Ihvw0rHqa8nm309uwm87aDbizi+hrrad4+mWkjxhD646VtGx+A5sjPlw30HmI8Tf/E6kFI9n253/Hao8PvXVl5jP6ynsHbaNl8+uEejuwOd342xtIL6mmePpCYuEAB1f8BWdaFkYshq9lP5WX3InNlUr9W0/RtnMV0z75H3iyjz4xs/6tp2jb9Q4VF99B7tjZg7YV6GqhY886yubedIb2WHJFAhFq/raX8bdUJ6oRdNZ2sfona7n4OxeR2j8Szdvk5dCGFqx2K20725nxuWnEwjHW/GI97qz4sPZQT5iR88sG3V9zaEMLnbVdONOc9LX24c52M+6msRgxg91P12D32AEL3XXdjLqyiqzKLGqe28OBV+s4/6uzyKvOTayrfkUDu56uYfxt1ZRfVDbo84R9Yfa/Wse4RWOx2s7+/oKCR+R9qFv2J4pnXIHN4SLs66LswltwpmbiPbSfA288ypS7vo0nu5hxN/4jzrRsTMOg5oX/IrUgPuIpq2IyZRcsOuH6oyE/rTtWMuWubwMQCfg4tGEJAK073sKTM4LiGZcDcHjbcqKhIJ6cEeRPuAB/eyNNaxcz+sr4Rf9osA9fy35c6TnHDR2AlNwSetw7ad2+koJJ807Xbho2O57YReXCikToRPwRmtY2484++ohvI2aw66kaZn15JharheKZRdjddixWC1PunkhufzBseWQbI2YXD9rG7r/tYcG/z8NisWDEDHb+ZTcAbdvjD9cbe/0YIP747UhfBJvDSsWCcjr3dFK7eF8ieIyYQXtNB3a37bihA+DOdJM7JofaxfsYe8OY07SXhs/ZH50iSRYLB+ht3J0IkcIpC3CmZsYnmgY2hxuA9BGjcaZlA9C1fxPZVdMS6zgSDvWrnsHbXDtoG1abg7Cvm9YdbxELB3F40ii/6FYAHCmZtO16G29zLaZpUjj5YtKLj1YnKJg0j56DOxJPFz28bTkFk+Ynpkf8vex9+dc0rX2RvS8/TKCrBYDsyqkc3rrsdO2mYRMNRunc20lm2dGbZ2ue28vY60YPmK+7rgesUPfGQfYu3kdvoxeb04YzzZkInVBvCCzgynDxbla7lQOv1RH2hbHarEz6eLyqtCvDyeFtbbTvaseIGRRMyqdoemFiuYLJ+UT6InTu6wKgeV0LI847GmyxcIwtf9jG3pf2sfWP2+nc2wlA3oRcGt5u4lyg4BE5RcGeNuzutONOa9nyBuXzbh/0fvvu1eSNm5P4u2T2tZTMvpbS869j36u/I9TbMWB+q93B+Ju/SnfdVjb+9mvsevan9LUeBCCvejZFUy+l7s0/s+Hhr1C3/AlikdCAZYumL6Rp7WJi4WA8uI4EI9DbtJdYKEDR9Csou2BRIigdqZkEOpvf/475kPC3+XGkHq1GUP9WAwWT8weFR7AzSPf+HsouKmX0NVXUvVlPV38YJJZd2Uj5hcd/dPmsL83Ad6iP5d9+i3d+vIb2mvh3mFWZxfhbqql9eT+v37+MrX/cPugm1VFXV1G7eB+mYdLb2EtG2dHHb/ta+vA2e6lYUM64RWNx9rfb5rBhRA3CfUO/4fXDSsEjcoqMaLT/kdED1a96hrzq2QN6HxAvEupMz02UwgFIK4rPY7U7Sckrxdeyf9D6UvPLGXvtF5n5+Z+TXTWNmuf/X2JawaR5TL7zW0y56zsEOptpXvfSgGWLpl5C575NNK5+jsLJ8wdMy66aSvqI0ez86w+pf+tpLLb4Z7FYbRixs798ixE14o+M7te1r5veRi+1S/YTCUSpX9FA1/5u7B47qUWp8dNrFgvZVZl07utOLGcaJj0He8gelX3c7aTkpTDlnkksfPASKhdWsOG/NxENxfdf0bRC5nx1Nhf/+zxM02TXU7sHLFs8owh/e4Cav+0Z0BsCyCzPoOzCUtb/ciNb/rANq/XoZ7HYLBiRs3+ggYJH5BS50rOJhgY+N6V+1dNklIwhu2oanfs2DZjWsmUZRVMvTfwd6GymdfuKxN/B7jbcWQWDtrNv6e8AsNrs5I6dnQiF1h0r6W3aA4AzLYvMsvGDAsPm9FA45RIiAR+ujLwB0wIdTeSNv5DJd34LV0YObTtXARAL+ROnBs9m7mw3EX8k8ffUv5vM6KuqGH1VFQ6PnfL5ZWRXZZFVkUnEF06MfAt0BkktOPq0z5bNhymcOvh7OWLro9sBsFgtFEzOx2qPB0RHTQeN/SV5nKlOCiblY0QHhoXFamH01VV01/WQ865g87f7yRmdzdyvn0/B5HwOvF4HxEfpGWHjuKf9zjYaXCByipxp2dhdKUQCPhyeNA5tXErr9hV4m/bQtPZFIgEvOaOmAxALB4n4ewcEi83poXPfZsK+LiIBH7ljZyWuFx0rFg5St/wJbA4Xga5DjLr8UwB4sotp3rCEnoM7MGIRQr3tVF5yN9FgH2073yYW8pNZNp6yuTcC8YEKbTvfJuTtpGPvhkSbPTkjCPu6KZy8AIjXlDsXHp/gznLj8DgI94VxHnPKbe/ifUQCUQ6+WU/FZRWk5HoYe8MYdjy5C0eqA1ema0DQNK87xLRPTznhdqwOa3zZFDv+9gDjbq7G7rLjznZTt6yevtY+TMOk73Af424ai2ma1L1ZT299L72NXkrnllA6twSAujfriQZjNLzdSHZVFntf2kdmWQZ9rX7K58VP9fmafeSOy8FiPfsL8alkjki/UymZ42s5QNuut6m85K4z3KrkiIUD7H/tESov/QR2d+qQl/uwlszpruuhaU0zE+8YP9xNOS2MmMG2x3Yw5rrRpPQ/vuHD5lRK5qjHI/I+pBVVYvekEQsHsTndw92cDywS8FG18JPnxGcByKrITDxywO4++w9zkb4I426uxpU+9Ec4fJid/d+IyDBxZ+YPdxNOm3Ppsxxx5Jk754Jz4brOsTS4QEREkkrBIyIiSfW+BheIiIi821AHF5xy8IiIiHwQOtUmIiJJpeAREZGkUvCIiEhSKXhERCSpFDwiIpJUCh4REUkqBY+IiCSVgkdERJJKwSMiIkml4BERkaT6/60+HRiltPzfAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 504x311.496 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "width = 7\n",
    "height = width / 1.618    # golden ratio\n",
    "fig, ax = plt.subplots(figsize=(width, height))\n",
    "\n",
    "squares = {'identical': len(psms_match),\n",
    "           'conflicting\\nsimilar': num_high_sim,\n",
    "           'conflicting\\ndifferent': num_low_sim,\n",
    "           'unique ANN-SoLo': len(psms_unique_annsolo),\n",
    "           'unique iPRG2012 consensus': len(psms_unique_consensus)}\n",
    "squares = {f'{key}\\n({value} SSMs)': value\n",
    "           for (key, value) in squares.items()}\n",
    "\n",
    "colors = sns.color_palette('Set1', len(squares))\n",
    "\n",
    "squarify.plot(sizes=squares.values(), color=colors, label=squares.keys(),\n",
    "              ax=ax, alpha=0.8, linewidth=1, edgecolor='black')\n",
    "\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "\n",
    "plt.savefig('iprg2012_consensus_treemap.pdf', dpi=300, bbox_inches='tight')\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "peptides_library = set()\n",
    "filename_pepidx = '../../data/interim/iprg2012/human_yeast_targetdecoy.pepidx'\n",
    "with open(filename_pepidx) as f_in:\n",
    "    for line in f_in:\n",
    "        if not line.startswith('#'):\n",
    "            peptide, info, _ = line.split()\n",
    "            charge = info.split('|')[0]\n",
    "            peptides_library.add((peptide, float(charge)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# conflicting iPRG consensus identification ions not found IN library: 58 (85 spectra)\n",
      "# conflicting iPRG consensus identification ions not found NOT IN library: 466 (657 spectra)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "psms_notfound = psms_different  # pd.concat([psms_different, psms_unique_consensus])\n",
    "peptides_notfound = set(zip(psms_notfound['sequence_consensus'],\n",
    "                            psms_notfound['charge_consensus']))\n",
    "peptides_notfound_library = peptides_notfound & peptides_library\n",
    "ions_consensus = psms['sequence_consensus'] + psms['charge_consensus'].map(str)\n",
    "psms_notfound_library = (psms[ions_consensus.isin(\n",
    "    [f'{seq}{charge}' for seq, charge in peptides_notfound_library])]\n",
    "                         .merge(psms_notfound))\n",
    "peptides_notfound_notlibrary = peptides_notfound - peptides_library\n",
    "psms_notfound_notlibrary = (psms[ions_consensus.isin(\n",
    "    [f'{seq}{charge}' for seq, charge in peptides_notfound_notlibrary])]\n",
    "                            .merge(psms_notfound))\n",
    "print(f'# conflicting iPRG consensus identification ions not found IN library: '\n",
    "      f'{len(peptides_notfound_library)} '\n",
    "      f'({len(psms_notfound_library)} spectra)\\n'\n",
    "      f'# conflicting iPRG consensus identification ions not found NOT IN library: '\n",
    "      f'{len(peptides_notfound_notlibrary)} '\n",
    "      f'({len(psms_notfound_notlibrary)} spectra)\\n')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
